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Recovery of the internal orbital structure of galaxies 
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ABSTRACT 

We construct axisymmetric and triaxial galaxy models with a phase-space distribution func- 
tion that depends on linear combinations of the three exact integrals of motion for a separa- 
ble potential. These Abel models, first introduced by Dejonghe & Laurent and subsequently 
extended by Mathieu & Dejonghe, are the axisymmetric and triaxial generalisations of the 
well-known spherical Osipkov-Merritt models. We show that the density and higher order 
velocity moments, as well as the line-of-sight velocity distribution (LOSVD) of these mod- 
els can be calculated efficiently and that they capture much of the rich internal dynamics 
of early-type galaxies. We build a triaxial and oblate axisymmetric galaxy model with pro- 
jected kinematics that mimic the two-dimensional kinematic observations that are obtained 
with integral-field spectrographs such as SAURON. We fit the simulated observations with ax- 
isymmetric and triaxial dynamical models constructed with our numerical implementation of 
Schwarzschild's orbit-superposition method. We find that Schwarzschild's method is able to 
recover the internal dynamics and three-integral distribution function of realistic models of 
early-type galaxies. 

Key words: stellar dynamics - celestial mechanics - galaxies: elliptical and lenticular, cD - 
galaxies: kinematics and dynamics - galaxies: structure 



!_< \ 1 INTRODUCTION 



The equilibrium state of a collisionless stellar system such as an 
elliptical or lenticular galaxy is completely described by its distri- 
bution function (DF) in the six-dimensional phase space of posi- 
tions and velocities. The recovery of the DF from observations is 
difficult, as for galaxies other than our own, we can usually only 
measure the projected surface brightness and the line-of-sight ve- 
locity distribution (LOSVD) of the integrated light as a function 
of position on the plane of the sky. Moreover, we generally do not 
know the intrinsic shape of the galaxy, nor the viewing direction, or 
the contribution to the gravitational potential provided by a super 
massive central black hole and/or an extended halo of dark matter. 
By Jeans (1915) theorem, the DF is a function of the isolating inte- 
grals of motion admitted by the potential, but it is not evident how 
to take advantage of this property other than for the limiting case of 
spherical systems. Orbits in axisymmetric geometry have two exact 
integrals of motion, the energy E and the angular momentum com- 
ponent L z parallel to the symmetry z-axis, but the third effective 
or non-classical integral Iz obeyed by all regular orbits is gener- 
ally not known in closed form. In stationary triaxial geometry E 
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is conserved, but regular orbits now have two additional effective 
integrals of motion, I2 and I3, which are not known explicitly. 

Schwarzschild (1979, 1982) devised a numerical method 
which sidesteps our ignorance about the non-classical integrals of 
motion. It allows for an arbitrary gravitational potential, which may 
include contributions from dark components, integrates the equa- 
tions of motion for a representative library of orbits, computes 
the density distribution of each orbit, and then determines the or- 
bital weights such that the combined orbital densities reproduce 
the density of the system. The best-fitting orbital weights represent 
the DF (cf. Vandervoort 1984). Pfenniger (1984) and Richstone 
& Tremaine (1984) included kinematic moments in this method, 
and Rix et al. (1997) showed how to include observed LOSVDs. 
A number of groups have developed independent numerical imple- 
mentations of Schwarzschild's method for axisymmetric geometry 
which fit the projected surface brightness and line-of-sight velocity 
distributions of early-type galaxies in detail (van der Marel et al. 
1998; Cretton et al. 1999; Gebhardt et al. 2000, Valluri, Merritt & 
Emsellem 2004; Thomas et al. 2004; Cappellari et al. 2006). Ap- 
plications include the determination of central black hole masses 
(see also van der Marel et al. 1997; Cretton & van den Bosch 
1999; Verolme et al. 2002; Cappellari et al. 2002; Gebhardt et al. 
2003; Copin, Cretton & Emsellem 2004), accurate global dynam- 
ical mass-to-light ratios (Cappellari et al. 2006), as well as dark 
matter profiles as a function of radius (Cretton, Rix & de Zeeuw 
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2000; Thomas et al. 2005), and recovery of the DF (Krajnovic et 
al. 2005). Van de Ven et al. (2006) and van den Bosch et al. (2006) 
included proper motion measurements in order to model nearby 
globular clusters, and determine their distance, inclination as well 
as mass-to-light ratio as function of radius. Finally, Verolme et 
al. (2003) and the companion paper van den Bosch et al. (2007, 
hereafter vdB07) describe an extension to triaxial geometry that in- 
cludes all line-of-sight kinematics. 

Although Schwarzschild models have significantly increased 
our understanding of the dynamical structure and evolution of 
early-type galaxies, questions remain about the uniqueness and the 
accuracy with which they are able to recover the global parameters 
as well as the internal dynamics of these galaxies. Many tests have 
been done to establish how the axisymmetric code recovers known 
input models, but these generally have been limited to spherical ge- 
ometry or to an input axisymmetric DF that is a function of E and 
L z only (van der Marel et al. 1998; Cretton et al. 1999; Verolme 
& de Zeeuw 2002; Valluri et al. 2004; Cretton & Emsellem 2004; 
Thomas et al. 2004; Krajnovic et al. 2005). 

One could construct a numerical galaxy model with 
Schwarzschild' s method itself, compute the observables, and then 
use these as input for the code and determine how well it recovers 
the input model. This is useful, but does not provide a fully inde- 
pendent test of the software. An alternative is to consider the special 
family of models with gravitational potential of Stackel form, for 
which all three integrals of motion are exact and known explicitly. 
These separable potentials have a core rather than a central cusp, so 
the corresponding models cannot include a central black hole, and 
are inadequate for describing galactic nuclei. However, they can be 
constructed for a large range of axis ratios (Statler 1987), and their 
observed kinematic properties are as rich as those seen in the main 
body of early-type galaxies (Statler 1991, 1994; Arnold, de Zeeuw 
& Hunter 1994). 

A small number of analytic DFs have been constructed for 
triaxial separable models. The 'thin-orbit' models (Hunter & de 
Zeeuw 1992) have the maximum possible streaming motions, but 
their DF contains delta functions, and they are therefore not par- 
ticularly useful for a test of general-purpose numerical machinery. 
Dejonghe & Laurent (1991, hereafter DL91) constructed separa- 
ble triaxial models in which the DF depends on a single param- 
eter S — E + wI-2 + Ills, which is a linear combination of the 
three exact integrals E, I2 and J3 admitted by these potentials, and 
is quadratic in the velocity components. For a given radial den- 
sity profile, the DF follows by simple inversion of an Abel integral 
equation. These so-called Abel models have no net mean streaming 
motions, and are the axisymmetric and triaxial generalisations of 
the well-known spherical Osipkov-Merritt models (Osipkov 1979; 
Merritt 1985), for which the observables can be calculated easily 
(Carollo, de Zeeuw & van der Marel 1995). Mathieu & Dejonghe 
(1999, hereafter MD99) generalised the results of DL91 by includ- 
ing two families of DF components with net internal mean motions 
around the long and the short axis, respectively, and compared the 
resulting models with observations of Centaurus A. Although the 
Abel character of the non-rotating components is no longer con- 
served, the expressions for the velocity moments in these more gen- 
eral models can still be evaluated in a straightforward way. When 
the entire DF depends on the same single variable S the famous 
ellipsoidal hypothesis (e.g., Eddington 1915; Chandrasekhar 1940) 
applies, so that self-consistency is only possible in the spherical 
case (Eddington 1916; Camm 1941). This does not hold for Abel 
models with a DF that is a sum of components for which the vari- 
able 5* has different values of the parameters w and u. Such multi- 



component Abel models can provide (nearly) self-consistent mod- 
els with a large variety of shapes and dynamics. 

Here, we show that for Abel models, in addition to the velocity 
moments, the full LOSVD can be calculated in a simple way. Next, 
we construct axisymmetric and triaxial Abel models to test our nu- 
merical implementation of Schwarzschild's method. We assume a 
convenient form for the gravitational potential, and construct the 
DF that reproduces a realistic surface brightness distribution. We 
compute the LOSVDs of the models and derive two-dimensional 
maps of the resulting kinematics. We show that, despite the simple 
form of the DF, these models display the large variety of features 
observed in early-type galaxies with integral-field spectrographs 
such as SAURON (Emsellem et al. 2004). By fitting axisymmetric 
and triaxial three-integral Schwarzschild models to the simulated 
observables we find that Schwarzschild's method is able to recover 
the internal dynamics and three-integral DF of early- type galaxies. 
In this paper we fix the mass-to-light ratio and viewing direction 
to those of the Abel models, while in our companion paper vdB07 
we investigate how well these global parameters can be determined 
by Schwarzschild's method, along with a full description of our 
numerical implementation in triaxial geometry. 

This paper is organised as follows. In Section[2]we summarise 
the properties of the triaxial Abel models of DL91 and MD99 
and present the intrinsic velocity moments in a form which fa- 
cilitates their numerical implementation. We describe the conver- 
sion to observables in Section [3] including the computation of the 
LOSVD. In Section|4]we construct a specific triaxial galaxy model 
and in Section [5] we fit the simulated observables with our triax- 
ial Schwarzschild models to investigate how well the intrinsic mo- 
ments and three-integral DF are recovered. In Section [6] we con- 
sider Abel models in the axisymmetric limit and construct a three- 
integral oblate galaxy model to test our axisymmetric implementa- 
tion of Schwarzschild's method. We summarise our conclusions in 
Section|7] In Appendix[A] we describe the simpler Abel models for 
the elliptic disc, large distance and spherical limit, and link them to 
the classical Osipkov-Merritt solutions for spheres. Readers who 
are mainly interested in the tests of the Schwarzschild method may 
skip Sectionsl2l-l4land[6~T1-[63l 



2 TRIAXIAL ABEL MODELS 

The triaxial Abel models introduced by DL91 have gravitational 
potentials of Stackel form, for which the equations of motion sep- 
arate in confocal ellipsoidal coordinates. We briefly describe these 
potentials, and refer for further details to de Zeeuw (1985a). We 
then make a specific choice for the DF, for which the velocity mo- 
ments simplify. 



2.1 Stackel potentials 

We define confocal ellipsoidal coordinates (A, /1, v) as the three 
roots for r of 
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r + a t + j3 r + 7 



(2.1) 



with (x, y, z) the usual Cartesian coordinates, and with constants 
a, /3 and 7 such that —7 < v < —f3 < n < —a < A. From the 
inverse relations 



2 _ (A + a)(p + a)(y + a) 
(a -/?)(<* -7) 



(2.2) 
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and similarly for y 2 and z 2 by cyclic permutation of a — *■ /3 — > 
7 — *• a, it follows that a combination (A, /i, i/) generally corre- 
sponds to eight different points (±x, ±y, ±2). In these coordinates, 
the Stackel potentials have the following form (Weinacht 1924) 



V s (X,p:,v) 



U(X) 



+ 



(A-/i)(A-i/) <ji-v){p-\) 

U(y) 



+ 



(v- X)(v- p)' 1 



(2.3) 



where U(t) is an arbitrary smooth function (r = X,p,v). The 
right-hand side of eq. d2.3t can be recognised as the second order 
divided difference of U (r). Henceforth, we denote it with the cus- 
tomary expression U[X, p, v], which is symmetric in its arguments 
(see Hunter & de Zeeuw 1992, eqs 2.1-2.3, 2.13 and 2.14). Addi- 
tion of a linear function of r to U(t) does not change Vs. 

The density ps that corresponds to Vs can be found from Pois- 
son's equation 



47rGp s (A, p, v) = V 2 V S (A, p, v), 



(2.4) 



or alternatively by application of Kuzmin's (1973) formula (see de 
Zeeuw 1985b). This formula shows that, once we have chosen the 
confocal coordinate system and the density along the short axis, 
the mass model is fixed everywhere by the requirement of separa- 
bilitjQ For centrally concentrated mass models, Vs has the s-axis 
as long-axis and the z-axis as short-axis. In most cases this is also 
true for the associated density (de Zeeuw, Peletier & Franx 1986). 

2.2 Orbital structure 

The Hamilton- Jacobi equation separates in (A, p, v) for the poten- 
tials J2.3I >, so that every orbit has three exact integrals of motion 
(cf. de Zeeuw & Lynden-Bell 1985) 

E = 1 {v't+v 2 y + vl) +U[X,fjt,u], 
h = \TL 2 y + \L 2 z + \{ a -l3)v 2 x 

+ (a- f3)x 2 U[X,p,v,-a], (2.5) 
h = \L\ + ±(l-T) J L 2 , + §( 7 -/?)u 2 

+ (l-P)z 2 U[X,v,v, -7], 



where v x 



and v z are the velocity components in the Carte- 



sian coordinate system, and L x 



the component 



of the angular momentum vector parallel to the x-axis. The other 
two components, L y and L z , follow by cyclic permutation of 
x — > y — > 2 — > x and v x — » v y — > v x — > v x . Furthermore, T 
is a triaxiality parameter defined as 



r=GS-a)/( 7 -a), 



(2.6) 



and U[X, p, v, a] is the third-order divided difference of U (r). All 
models for which U"'(t) > have a similar orbital structure and 
support four families of regular orbits: boxes with no net rotation, 
inner and outer long-axis tubes with net rotation around the a>axis, 
and short-axis tubes with net rotation around the z-axis (Kuzmin 
1973; de Zeeuw 1985a; Hunter & de Zeeuw 1992). 



A third method for the calculation of the density is to use 4nGps = 
H[X, A, fj,. fj,, v, u], where the fifth-order divided difference is of the func- 
tion H(t) = 4ci(t)U'(t) - 2o'(r)E/(r) with o(t) = (r + a)(r + 
0)( T + 7) and U{t) defines the potential as in eq. (375). This result was ob- 
tained by Hunter in 1989 (priv. comm.) and by Mathieu & Dejonghe (1996). 
Similar expressions exist for the related families of potential-density pairs 
introduced in de Zeeuw & Pfenniger (1988). 



According to Jeans (1915) theorem the phase-space distribu- 
tion function (DF) is a function f(E, I2, 13) of the isolating inte- 
grals of motion (cf. Lynden-Bell 1962; Binney 1982). The velocity 
moments of the DF are defined as 



Hlmn(X, fJ,, V) 



III I m n 

JJJ V ^ Vv 



/(B,Ja,l3)d«Ad« M d«„, (2.7) 



where /, m and n are non-negative integers, and v\, « M and v v are 
the velocity components in the confocal ellipsoidal coordinate sys- 
tem. Many of the velocity moments vanish due to the symmetry of 
the orbits in these coordinates. The zeroth-order velocity moment 
is the mass density that corresponds to the DF 



(2.8) 



/0*(A, n, v) = |Uooo(A, fi, v). 

In self-consistent models, p + must equal ps given in eq. d2.4| l, the 
mass density that is related to the potential Vs by Poisson's equa- 
tion. 



2.3 Abel distribution function 

Following DL91, we choose the DF to be a function of the three 
integrals of motion E, I2 and 73 as given in eq. i2.5\ through one 
variable 

f(E,I 2 ,I 3 )= f(S), with S= -E + wh+uh, (2.9) 

and w and u are two parameter^ This choice for the DF is 
equivalent to the celebrated ellipsoidal hypothesis (e.g., Edding- 
ton 1915; Chandrasekhar 1940). Self-consistency is only possible 
in the spherical case (Eddington 1916; Camm 1941). On the other 
hand, these DFs can produce realistic (luminous) mass densities p*, 
which differ from the (total) mass density ps, as in galaxies with 
dark matter (see also § 12.41 below when we combine DFs of the 
form 12.91 with different values for w and u.) 

DL91 and MD99 divided the DF into three types of compo- 
nents. The non-rotating (NR) type is made of box orbits and tube 
orbits with both senses of rotation populated equally. The two rotat- 
ing types, LR and SR, consist of tube orbits, and have net rotation 
around either the long axis or the short axis. 

2.3.1 Velocity moments 



Due to the choice \2.9l of the DF, the general expression l !2.7t for 
the velocity moments can be simplified, as shown by DL91 for 
the non-rotating components and by MD99 for the rotating compo- 
nents. We recast their expressions into a different form to facilitate 
the numerical implementation. The resulting velocity moments are 
given by 



2l + m+n+3 
pimn(X, p,u) = , rrm + 1 rrn + 1 



xj T lmn [S top (X,pt,u)-S] (l+m+n+1)/2 f(S)dS, (2.10) 

•Smin 

and set to zero at positions for which S ma x < Smin- The terms 
H\i,v, H v \ and Hx^ in the square root in front of the integral are 
defined as 

H„ T = 1 + ( " ± a)(T + a) w+ {a + 7)(T + 7) u, (2.11) 

7 — ct a — 7 

2 In contrast with DL91 and MD99, we choose Vs < and E < 0, 
consistent with e.g. de Zeeuw (1985a). 




Figure 1. The limiting value Sy lm of the variable S = — E + WI2 + UI3 
as function of the parameters to and u. The physical region is bounded by 
the relations j2.12K indicated by the thick solid lines. The dashed curves 
divide this region into three parts, each with a different expression for Sii m . 
The relations for the separatrices L\ and L2 are given in eq. 12.13) . 



with a, t — X,fi,i>. Orbits are confined to the region of space for 
which all three terms are non-negative. In general, this condition 
will not be satisfied for all points, so that the Abel components 
have finite extent. From the requirement that at least the origin 
(A, fj,, v) = (—a, — P, —7) should be included, we find the fol- 
lowing limits on w and u 



w > 



and u < 



y-0 



(2.12) 



The factor Ti mn in the integrand as well as the upper limit S max 
of the integral are different for each of the three Abel component 
types NR, LR and SR, and are discussed in §§ 12. 3. 2j|2. 3.41 below. 
The lower limit of the integral S m ln has to be at least as large as the 
smallest value possible for the variable S. This limiting value Sum 
depends on the choice of the DF parameters w and u in d2.9t . as is 
shown in Fig.[T](cf. Fig. 7 of DL91). The boundaries follow from 
( |2.12t and the separatrices L 1 and L2 are given by 



Li 



u[-ot, z -7,-7] 

09 -a) 
u 

1 — (7 — a) u 



(2.13) 



At a given position (A, fj,, v), orbits with different values of the in- 
tegrals of motion E, I2 and I3, and hence different values of S, 
can contribute to the integral ( 12. 1 0> . The restriction to bound or- 
bits (E < 0) together with the requirement that v\, w„ and v„ all 
three have to be non-negative determines the part of the integral 
space that is accessible by orbits that go through (A, fj,, v). An ex- 
ample of the resulting tetrahedron in the [E, I2, /3)-space is shown 
in Fig. |2](cf. Fig. 1 of MD99). The largest possible value of S is 



Figure 2. The tetrahedron shows all accessible points in integral space 
(E, I2, 1 3) for a given position (A, fj,, v). The tetrahedron is bounded by 
the planes for which v\ = 0, = 0, v^, = and E = 0, respectively. 



The two shaded planes, which are given by vi 



at A : 



I-' - 



and = = at fi = v = — (3, divide the tetrahedron into the 
parts corresponding to the four general orbit families in a triaxial separable 
potential: box (B) orbits, inner (I) and outer (O) long-axis tube orbits and 
short-axis (S) tube orbits. 



given by the top of this tetrahedron 

1 Sto P (A, J u, v) = — U[X,/x, v] 

(A + a)((ji + a)(y + a) 



7 — a 
(A + 7)^ + 7)^ + 7) 

a — 7 



U[X,n,v,-a] 
U[X,H,u,-i\, (2.14) 



which is thus a function of the position (A,/i, v). At the origin 
Stop (—a, — j8, —7) = U[— a, —/3, —7], which is the central value 
of the potential Vs. In what follows, we normalise Vs by setting 

U\-a, -0, -7] = -1, so that < S top < 1. 



2.3.2 Non-rotating components (NR) 

Since the non-rotating component type can exist everywhere in the 
accessible integral space (the tetrahedron in Fig.[2}, we simply have 
that S m ax = Stop (A, n, v). Spatially the NR components are thus 
bounded by the surface Stop (A, p, v) = S m in. 

The factor Ti mn follows from the cross section of the S-plane 
within the tetrahedron and can be written in compact form as 



Tr^ n = B(^,m±l,n±^), (2.15) 

where B is the beta function of three variable^]- Since is in- 
dependent of S it can be taken out of the integral (cf. eq. [3.10] of 

3 The beta function of k variables in terms of the complete gamma function 
T is defined as B(/3i ,...,&) = T(/3i) ■ ■ ■ V{p k )/V{fli +•••+&). 
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DL91), which then becomes of Abel form. Unfortunately, the inver- 
sion of eq. J2. lOt for any chosen moment p; mn (A, p, v), including 
the case I — m — n — 0, is generally impossible, as the left-hand 
side is a function of three variables, while the DF depends on only 
one variable, S. The density p* specified along any given curve will 
define a different f(S). A case of particular interest is to choose the 
density along the short axis to be p*(0, 0, z) — ps(0, 0, z). This 
defines a unique f(S), and hence gives p* everywhere. Kuzmin's 
formula applied to ps(0, 0, z) similarly defines the density ps ev- 
erywhere. For single Abel DF components these will not be the 
same, except in the spherical limit (see Appendix lA3t . 

Since the orbits have no net rotation, the velocity moments 
l-ifmn are on ly non-zero when I, m and n are all three even, and 
vanish in all other cases. 



2.3.3 Long-axis rotating components (LR) 

The long-axis rotating component type only exists in the part of 
the integral space that is accessible by the (inner and outer) long- 
axis tube orbits. Within the tetrahedron for all orbits this is the 
region for which v'l > at v = —/3. It follows that S ma x = 
St op (A,p, — 0) < S top (A,p, v). 

The term T; mn follows from the cross section of the S-plane 
within the tetrahedron and with the above boundary plane v% = 
at v — —/3. Without any further constraint this results in zero net 
rotation, because each orbit with positive rotation around the long 
axis with v v > 0, is balanced by an orbit with opposite direction of 
rotation with v v < 0. Therefore, we restrict to orbits with v v > 0, 
resulting in maximum streaming around the long axis for each LR 
component. This reduces the accessible integral space, and thus 
also the term Timn, by a factor of two, so that the latter becomes 



! (_2)a+™)/y^+ 1 6™+ 1 M& 1 



(s + l)(s-l)...(s + l-(J + m))' 
with s = I + m + n, the parameters ao and bo defined as 

(A + /3)# M „ [Stop(A, p, —0) — S] 



(2.16) 



£1(1 



bo = 



(A - v) i? M (-/3) [Stop (A, p, v) - S] ' 



(2.17) 



(p - v) H(-@)x [S top (A, p, v) - S] ' 
which for S < S max = Stop (A, P, —0) are non-negative, and 



M(s 



I 

2 ' ' 



; ao, bo, ■ 



M(s ,f i;6 ,o 0) f) 



a < b , 
ao > b Q . 



(2.18) 



The function A4(s,i, j; a,b, <p) is defined in Appendix iBl where 
we evaluate it in terms of elementary functions (odd s) and elliptic 
integrals (even s). 

The LR components have maximum streaming around the 
long axis, but the motion parallel to the intermediate axis and short 
axis cancels. As a result, the velocity moments p^ n vanish when / 
or m are odc0. Multiplying pf^ n with (—1)" results in maximum 
streaming in the opposite direction. By choosing different weights 
for both senses of rotation, we can control the direction and the 
amount of long-axis streaming motion for each LR component. 



4 Since I + m is even, the factor (— 2)('+ m )/ 2 in eq. 12.161 is always real. 



2.3.4 Short-axis rotating components (SR) 

The short-axis component type reaches the part of integral space 
accessible by the short-axis tube orbits. Within the tetrahedron for 
all orbits this is the region for which v'^ > both at p = —f3 and 
p = —a (Fig.[2](. The latter requirement is equivalent to I2 > 0. In 
this case, Sm^ = Stop(A, — a, v ) < 5 , top (A, p, v). 

The form of the term T; mn depends on the cross section of 
the S-plane within the tetrahedron and with the above two bound- 
ary planes. In case each SR component has maximum streaming 
around the short axis (« M > 0), it is given by 



2 (-2) 



(i+n)/2 



(2.19) 



+ l)0-l)...(s + l- (l + n)) ' 

The parameters ai and Cr follow from ao and feo defined in ( 12. 17t 
by interchanging p «-»• v, and in turn 02 and C2 follow from a\ 
and Ci by interchanging a <-» (i. For the terms A^f R we have two 
possibilities, I and II, 



Mf K = 



AC 



' M(s, |, f ;ai,ci,9x), ai < ci, 

M(s, |,|;ci,ai, f) 

-M(s, |-,|;ci,ai, f-ftt), ai > ci, 

'M(s, |,f 5011,011,1) 



(2.20) 



-M(s, 



;au,cn,6i{), on < cn, (2.21) 



K M(s, ~, |;cn,an, | — On), an > cn, 
where A4 is given in Appendix[B] and 8\ and 9n follow from 

cn (ai — an) 



tan 9j = 



an (cn — ci) 



, 2, ci (a n - 01) 
tan &n = — 7 f- 

ai (ci - cn) 



(2.22) 



For the assignment of the labels I and II, we discriminate between 
four cases 



ai < a 2 , Ci > C2 

ai > a 2 , ci < C2 

ai < a 2 , Ci < C2 

ai > a 2 , ci > C2 



II 
II 



2, 
1, 

tt/2, (7 2 si 
t/2, Cf 



0, 
0. 



(2.23) 



The SR components have maximum streaming around the short 
axis, so that the velocity moments pf„ n vanish when / or n are 
odd. Multiplying pf^n w i m (~l) m results in SR components with 
maximum streaming around the short axis in the opposite direction. 



2.4 Combination of multiple DF components 

Until now, we have chosen the Abel DF to be a function of a sin- 
gle variable S = —E + 111I2 + 11I3, and we have separated it in 
three component types, NR, LR and SR, but we have not made 
any assumption about the form of the DF (apart from the obvious 
requirement that it has to be non-negative everywhere and that it 
decreases to zero at large radii). Following MD99, we choose the 
DF to be a linear combination of basis functions of the form 



fs(S) = 



S Smiii 



(2.24) 



1 Smin 

which, like the velocity moments l |2.10t , are non-vanishing as long 
as Slim < S m in < S < S m ax < Stop < 1. The exponent S is a 
(non-negative) constant. 

Once the Stackel potential l |2.3t is known by defining the func- 
tion U(t), we can use the above relations (§[23) together with the 
expressions in Appendix [B] to compute for a given basis function 
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fs(S) the velocity moments l !2.10t for the NR, SR and LR compo- 
nents in an efficient way, where at most the integral over S has to 
be evaluated numerically. For the NR components this integral can 
even be evaluated explicitly, resulting in 



NR I x \ 



[2{S„ 



5*min)] 



i + m+n + 3 



S(iti,afi,a±l,tf+1), (2.25) 



where 5 max = Stop (A, fi, v) (cf. eq. |2.14| >. 

Each DF component and corresponding velocity moments 
thus depend on the choice of the DF parameters w, u and 5, the 
type of component, and for the rotating components (LR and SR), 
they also depend on the sense of rotation around the axis of sym- 
metry. By summing a series of DF basis functions over w, u and 
5, one might even expect to cover a large fraction of all physical 
DFs. Due to the different values of w and u, such a sum of DF 
components is no longer a function of the same, single variable 
S, so that the ellipsoidal hypothesis does not apply. Consequently, 
it becomes possible to construct (nearly) self-consistent dynamical 
models, with the (combined) luminous mass density p* equal (or 
close) to the mass density ps associated to the potential. 



3 OBSERVABLES 

We describe how the intrinsic velocity moments can be converted 
to projected velocity moments on the plane of the sky. Alterna- 
tively, these line-of-sight velocity moments follow as moments of 
the LOSVD, which we show can be calculated in a straightfor- 
ward way for Abel models. Parameterising the LOSVD as a Gauss- 
Hermite series, we obtain the observable quantities: the surface 
brightness, the mean line-of-sight velocity V, velocity dispersion 
a, and higher-order Gauss-Hermite moments /13, Ha, . . . 



3.1 From intrinsic to observer's coordinate system 

In order to calculate line-of-sight velocity moments, we introduce 
a new Cartesian coordinate system (x", y" , z"), with x" and y" in 
the plane of the sky and the z"-axis along the line-of-sight. Choos- 
ing the x"-axis in the (x,y)-plane of the intrinsic coordinate sys- 
tem (cf. de Zeeuw & Franx 1989 and their Fig. 2), the transforma- 
tion between both coordinate systems is known once two viewing 
angles, the polar angle 1? and azimuthal angle tp, are specified. The 
intrinsic z-axis projects onto the j/"-axis, which for an axisymmet- 
ric galaxy model aligns with the short axis of the projected surface 
density E. However, for a triaxial galaxy model the j/'-axis in gen- 
eral lies at an angle tp with respect to the short axis of E. This 
misalignment ip can be expressed in terms of the viewing angles $ 
and tp and the triaxiality parameter T (defined in eq. l2.6t as follows 
(cf. eq. B9 of Franx 1988) 



tan2t/i 



T sin 2<p cos ■d 



sin 2 § — T (cos 2 tp — sin 2 ip cos 2 $) 



(3.1) 



with sin 2ip sin 2p cos 1? < and —n/2 < ip < n/2. A ro- 
tation over ip transforms the coordinate system (x" ,y" , z") to 
(x',y', z'), with the a;'-axis and y'-axis aligned with respectively 
the major and minor axis of E, whereas z' = z" is along the line- 
of-sight. 

The expressions in § I2.3I involve the velocity components in 



the confocal ellipsoidal coordinate system The conver- 

sion to line-of-sight quantities can be done by four successive ma- 
trix transformations. First, we obtain the velocity components in the 
first octant of the intrinsic Cartesian coordinate system (x, y, z) via 
Q, of which the first element is given by (cf. DL91) 



Qi 



■ sign(A + a) 



(jx + a){v + a)(X + f3)(\ + 7) 



(3.2) 



and the other elements follow horizontally by cyclic permutation 
of A— * v — »A and vertically by cyclic permutation of 
q — ► /3 — ► 7 — > The second matrix uses the symmetries of the 
orbits to compute the appropriate signs of the intrinsic Cartesian 
velocities in the other octants. The result depends on whether or 
not the orbit has a definite sense of rotation in one of the confocal 
coordinates. For the three types of Abel components this results in 
the following matrices 



NR : S = diag[sgn(a;),sgn(y),sgn(«)] 
LR : S = diag[sgn(a;y2:),sgn(z),sgn(y)] 
SR : S = diag[sgn(y),sgn(a;),sgn(a;yz)] 



(3.3) 



Finally, the conversion from the intrinsic to the observer's Carte- 
sian velocities involves the same projection and rotation as for the 
coordinates. We represent these two coordinate transformations re- 
spectively by the projection matrix 



(— sin p cos tp i 

— cos 1? cos tp — cos $ sin p sin 1} 
sin & cos p sin $ sin tp cos 1} J 

and the rotation matrix 

(cos tp — sin ip 0\ 
sin?/) cos'0 I . 
l) 

In this way, we arrive at the following relation 



with M = RPSQ, 











|=M| 











(3.4) 



(3.5) 



(3.6) 



where the full transformation matrix M is thus a function of 

(A, fi, v), the constants (a, j3, 7) and the viewing angles (§, tp, ip). 

3.2 Line-of-sight velocity moments 

We can now write each velocity moment in the observer's Carte- 
sian coordinate system (x',y', z') as a linear combination of the 
velocity moments in the confocal ellipsoidal coordinate system 



fiijk(x',y',z') = 2j Cl,m,n{s) fJ,lmn(X, fJ,, v), 

l.m.n 



(3.7) 



with s — i + j + k = l + m + n. The coefficients c; im , n (s) are 
products of elements of the transformation matrix M in eq. i3.6l . 
They can be obtained with the following recursive algorithm 

!ci,o,o(s) ci-i >m ,n(s - 1), if i > 0, 
Co,i,o(s) co, m -i,n(s - 1), if m > 0, (3.8) 
co,o,i(s) co,o,n-i(s - 1), ifn>0, 

with the first order expressions given by 



ci,o,o(s) = M e „i, c ,i,o(s) = M e 



C0,0,l(s) = M e 



(3.9) 
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and Co, 0,0 = 1. The index e s is the sth element of the vector 
e = [3, .., 3, 2, .., 2, 1, .., 1], where the number of integers 3 (#3) 
is equal to the value of the velocity moment index k, and similarly 
#2 = j and #1 = i. 

The line-of-sight velocity moments now follow from (numer- 
ical) integration of ^loofc along the line-of-sight 



fik{x,y)= Vook(x ,y ,z )dz , (3.10) 

J — oo 

which are thus functions of position on the sky plane. 
3.3 Line-of-sight velocity distribution 

Using the definition of the intrinsic velocity moments of the DF 
(eq. 12. 7b and rearranging the sequence of integration, we rewrite 
eq. d3 - 1 01 > for the line-of-sight velocity moments as 

/oo 
v k z >C{x ,y' ,v z ,)dv z , , (3.11) 
- OO 

where we have introduced the LOSVD 

C(x',y',v,,) = Jjjf{E,h,h)&v x ,Av y ,&z'. (3.12) 

Although the integral over z' in general can only be evaluated nu- 
merically, we show that for the choice l |2.9b of the DF, the double 
integral over the velocities can be simplified significantly. 

Our analysis generalises the results for the well-known spher- 
ical Osipkov-Merritt models. We describe the spherical limit to- 
gether with the elliptic disc and large distance limit in Appendix A, 
while we present axisymmetric Abel models in § 16, 11 



3.3.1 Abel LOSVD 

Substituting the expressions l |2,5b for the integrals of motion in S = 

—E + WI2+ uLi, we obtain 

S = 5 top (A, fi,u)-± {Hfiv v\ + H vX vl + vl) , (3.13) 

where the expression for Stop (A, v) is given by eq. J2.14b and 
the terms H^, H v \ and H\^ L are defined in eq. J2. lib . Defining 



X 2 



2 [Stop (A, -S] 



■2 



and similarly Y and Z by cyclic permutation of A — * fj, 
we can write the expression for S as 



X 2 + Y 2 + Z 2 



1. 



(3.14) 

/ -> A, 

(3.15) 



For a given position (A, fj,, v), each value of S thus defines the sur- 
face of the unit sphere in the variables (X, Y, Z). In these variables, 
we can write the integral of the DF over velocities, i.e., the stellar 
mass density, as 



P* = 



/// /(s) 



dv x i dVyi dv z 



2 [Stop — S] 



f(S) 



jjJdXdYdZ 

X 2 +Y 2 +Z 2 =l 



dS. (3.16) 



This is the same expression as for the zeroth-order velocity moment 
of the DF, /iooo, in eq. ( 12.10b . where 2 Tboo is equal to the integral 
between square brackets. 

The matrix M in eq. ( 13 .6b provides the conversion from the 
velocity components in the confocal ellipsoidal coordinate system, 



(v\, Vp, Vv), to those in the observer's Cartesian coordinate system, 



,v z i). Hence, for a given line-of-sight velocity v z i , we find 



ei X + e 2 Y + e 3 Z = v z , / g(S). 

The coefficients ei, e 2 and e$ are defined as 

het = \J H ly \H\^ A/31, 



he 2 = \[Ey^M^ M32, 
he 3 = ^jH^Hvx M 33 , 
and normalised with respect to h given by 

h 2 — H V \H\^ + Hx^H^v M32 + Hp, v H v \ M| 3 



(3.17) 



(3.18) 



(3.19) 



These coefficients are functions of the position (A, /i, v), the con- 
stants (a, P, 7) and the viewing angles (#, ip, tp) through the com- 
ponents of the matrix M, and also depend on the DF parameters w 
and u through the terms H\fj,, and H v \. It follows that 



g(S) = h 



2[S U 



S] 



(3.20) 



which is a function of the variable S. 

We thus find that each combination of values of S and v z i 
results in the cross section of the surface of the unit sphere in 
eq. d3.15b with the plane in eq. ( 13.17b . i.e., a circle, in the vari- 
ables (X, Y, Z). We rotate the latter coordinate system such that 
the normal vector {e\,ei,ez) of the plane of the circle coincides 
with the Z'-axis of the system given by 



f X\ I cos$ sin $ cos 6 sin $ sin 6 
Y I = I — sin $ cos <E> cos cos $ sin O 
,Zj \ -sine cos6 




(3.21) 



where the rotation angles $ and follow from 

tan$ = - J tan0 = ^fl±£l. (3 . 2 2) 

In these coordinates the circle is conveniently parameterised as 

(3.23) 



X' = Vl ~ Z' 2 cos£', Y' = s/l- Z' 2 sing, 



where Z' — v z i/g(S). We can now rewrite the integral between 
square brackets in eq. < I3 . 16b as 



// 



<9R <m 

d? A dZ> 



d^' dZ' = jj d£ dv z , , (3.24) 



where the vector R = (X, Y, Z) and A indicates the cross product. 
The integral over £' is the length of the part of the circle, A£', 
for which the corresponding integral space is accessible by orbits, 
and hence is in general a function of S and v z i and differs for the 
different types of Abel components as we show below. 
Inserting eq. l |3,24b in eq. l |3.16b , we obtain 



9{S) 



f(S) A£'(tv,S)chv as, 



-s(S) 

tin.) s «p(v) 



f(S) M'(v z ,,S) dS dv z 



(3.25) 



-9(S mi „) 



where after changing the order of integration in the last step, the 
upper limit of S is given by S up — min[G(iv), S max ], with 



G(v z i) — S t0 p(A, fi, v) — H^HuxHxn^^- 



(3.26) 
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Comparing the first line of eq. J3,16t with the second line of 
eq. d3.25t . we see that the choice of the Abel DF, f(E, I 2 , 7 3 ) = 
/(S), indeed reduces the triple integration d3.12t for the LOSVD 
to a double integral: 



£(x',y',v z i - / - 



f(S)A£(v z ,,S)dSdz', (3.27) 



Smin 



and vanishes when \v z /\ exceeds the 'terminal velocity' v t = 
<?(S'min)- The expressions for h and S up follow from eqs d3. 19\ 
and \3.26h whereas S max and A£' are different for each of the 
three Abel component types and are considered next. 

3.3.2 Non-rotating components (NR) 



As for the intrinsic moments in § 12.3.21 we have for the non-rotating 
component type that Smax = Stop (A, /i, v), and, since the full in- 
tegral space is accessible, A£nr = 27r, independent of S and v z i. 

In the case of a basis function fs(S) as defined in eq. J2.24I ). 
the integral over S can be evaluated explicitly resulting in 



r NR _ 
'-S — 



2tt 



(5+l)(l-S n 



- [G(*v)-S min f +1 dz'. (3.28) 



3.3.3 Long-axis rotating components (LR) 

The integral space accessible by the (inner and outer) long-axis tube 
orbits is given by v 2 > at v — —(3, so that immediately S max = 
Stop (A, n, —0), whereas the calculation of A£l R is more complex. 

Since eq. d3.15t must also hold at v — —j3, we find that for 
LR components, within the unit sphere in the variables (X, Y, Z), 
the space is restricted to that within the elliptic cylinder given by 

X 2 Y 2 

— + — = 1, (3.29) 
a b 



(3.31) 



where ao and bo are defined in eq. d2.17t . In the rotated coordi- 
nate system (X',Y',Z') defined in eq. d3.2U , at height Z' = 
v z i I 'g(S), the elliptic cylinder results in an ellipse given by 

diX' 2 + d 2 Y' 2 + d 3 X'Y' + d 4 X' + d 5 Y' + d 6 = (3.30) 

for < £' < 27T and with coefficients 

d\ = (b el + a e 2 ), 

d 2 = e 2 (bo e 2 + a el), 

d 3 = 2eie 2 e 3 (6 - a ), 

d 4 = 2eie 2 (ef + el) 5 (b — ao) (v z >/g), 

d 5 = 2e 3 (e? + el)? (b a e{ + a el) (v z ,/g), 

da = (e? + el) [(b e( + a el) (v z >/g) 2 - a b ] . 

Because all the above relations only involve the squared values of 
X, Y and Z, they are independent of the sign of the corresponding 
velocities V\, u M and v v (cf. eq. 13.14b . which results in zero net 
rotation. For the LR components, to obtain net rotation around the 
long axis, we simply limit the range of v v values, e.g., requiring 
Z > 0, results in maximum streaming motion around the long 
axis. This restricts the space in (X r , Y' , Z'), at given height Z' — 
v z i I 'g(S), to one side of the line 



'{el + elY Y' + e z (v z ,/g) = 0. 



(3.32) 



around the long axis. By choosing different weights for both senses 
of rotation, we can control the direction and the amount of long- 
axis streaming motion. 

For given values of S and v z i, the integral space covered by 
the LR components is thus the part of the circle in eq. d3.23t that 
falls within the ellipse in eq. < |3 ,301 > and that is on the correct side 
of the line in eq. d3.32t (see also Fig. |3j- The length Ac;l R of this 
part thus ranges from zero to a maximum of 27r when the circle 
is completely inside the ellipse and on the correct side of the line. 
To compute this length, we determine the points where the circle 
(possibly) intersects the ellipse and the line. Substituting the circle 
parameterisation of eq. fl3.23t in the expression for the ellipse in 
eq. ( 13.30b . we find that the intersections with the ellipse are the 
(real) zero points of the following fourth order polynomial in u = 
tan(£'/2) 

(de - d 4 R' + dii?' 2 ) u 4 + 2R'(d 5 - d 3 R') u 3 

+ 2 [d 6 + (2d 2 - rfi)i?' 2 ] u + 2R'(d 5 + d 3 R') u 

+ (d e + d 4 R' + dii?' 2 ) = 0, (3.33) 



where we have introduced R' = y/1 — (v z i / g) 2 - The intersections 
with the line result in the following two solutions 



u± = 



(e 2 + e 2 )iR'±[l- 



e\ 



(v*>/9) 2 



ez (v z '/g) 



(3.34) 



for \v z i | < g(S)(l — e|) 2, otherwise the line is outside the circle. 

We thus (numerically) find up to six real zero points u% and 
corresponding angles ^ = 2 arctan(iii), sorted from low to high. 
For the set {— n, £' 2 , . . . , n}, we compute the lengths of the se- 
quential intervals on the circle for which the corresponding values 
fall within the ellipse and on the correct side of the line. This can be 
checked by inserting a value from the corresponding interval, e.g. 
the central value, in eq. J3.23I ) and substituting the resulting X' and 
Y' into eqs ( 13.301 ) and l !3.32t . If the left-hand side is negative (posi- 
tive), the interval is inside (outside) the ellipse, and (for Z > 0) on 
the wrong (correct) side of the line. Finally, the sum of the resulting 
interval lengths provides A£l R . 

3.3.4 Short-axis rotating components (SR) 

The short-axis tube orbits are restricted to the region of integral 
space for which v 2 . > both at /i = — j3 and fi — —a, and hence 
S m ax = Stop (A, —a, v). For the calculation of A£s R we have that 
the space within the unit sphere in (X, Y, Z) is now restricted to 
the part that falls within both the elliptic cylinders 



X 2 z 2 _ | 

ai ci 



and 



x 2 z 2 _ 1 

0,2 C2 



(3.35) 



As in § |2.3.4| for the intrinsic moments, a\ and c 1 follow from ao 
and 6o defined in ( 12. 17t by interchanging y, <-> v, and in turn a 2 
and c 2 follow from ai and ci by interchanging a <-» j3. 

Both elliptic cylinders result in ellipses in the Z'-plane, as in 
eq. < !3.30b for the LR components, but now with coefficients 



(3.36) 



The restriction to the opposite side of the line inverts the rotation 





= Ciel, 




d 2 


2 2 . / 2 . 2 \ 2 

= d e ± e 3 + ai (ex + e 2 ) , 




dz 


= 2cjeie 2 e3 




a*4 


= 2d eie 2 (e 2 + e\) 5 (v z < 


/<?), 


d 5 


= 2e 3 (el + el) z [a el - 


a-i (el + el)] (v z ,/g), 


do 


= (el + el) [(a el + ai e\ 


) ( v z'/g) 2 - am] . 




Figure 3. The part of integral space accessible by long-axis (left panel) and short-axis (right panel) tube orbits, for given values of the Abel DF variable S 
and the line-of-sight velocity v z i . This corresponds to the thick part of the circle, which is restricted to be within the dashed ellipse(s) and below [above] the 
dashed line for long-axis [short-axis] rotating components with maximum streaming (see text for details). The length of the thick part of the circle equals A£' 
in the expression )3.27t of the LOSVD. 



for i = 1 , 2 respectively. The zero points of the corresponding 
fourth order polynomials in eq. {333} are again the intersections 
with the circle in eq. J3.23t . 

Net rotation around the short axis follows by limiting the range 
of Vfi values, e.g., Y > yields maximum streaming, which re- 
stricts the accessible integral space to one side of the line 



-ex X' + e 2 e 3 Y 1 + e 2 (e? + el) ~- (v z ,/g) = 0. 
The two solutions of the intersection with the circle are 



u± = 



-e 2 e 3 ii'±(e? + e|)2 [l . 



el 



(vz'/g) 



21 2 



exR' + {e 2 1 +e 2 2 )2e 2 (iv/ff) 



(3.37) 



(3.38) 



for \v z / \ < g(S)(l — el) 2, otherwise the line is outside the circle. 

The combination of all the (real) zero points provides the (or- 
dered) set {— 7r, £x)£a! ■ ■ ■ i 71 "}' w ' m at m °st ten intersections ^ 
with the circle given in eq. <3.23t . We compute the lengths of the 
circle intervals for which the enclosed values fall within both el- 
lipses and on the correct side of the line. This means, for which the 
corresponding X' and Y' values substituted in eq. J3.30t result in a 
negative left-hand side for both pairs of a% and bi, and (for Y > 0) 
in a positive left-hand side of eq. 13.371 . Finally, A£g R is the sum 
of the resulting interval lengths. 



3.3.5 Other type of components 

When considering the LR type of components we make no dis- 
tinction between inner and outer long-axis tube orbits because they 
have similar dynamical properties. Similarly, the non-rotating box 
orbits are part of the NR type of components and are not considered 
separately. Nevertheless, if we are interested in the specific contri- 
bution of these orbit families to the LOSVD, this can be achieved 
by a straightforward extension of the above analysis. 

As can be seen from Fig. [2] the inner and outer long-axis tube 



orbits are separated by the plane 1% = 0, or equivalently the regions 
for which v\ > at A = —a and v 2 . > at /i = —a, respec- 
tively. This is in addition to the restriction v 2 > at v = — /3 for 
both long-axis tube orbits. For the inner long-axis tube orbits this 
implies that 5 max = -Stop (A, /x, — /?)■ The space within the unit 
sphere in (X, Y, Z) is now restricted to the part that falls within the 
intersection of the elliptic cylinders in eq. l !3.29t and 

Y 2 7 2 

— + — = 1, (3.39) 
03 c 3 



where 63 and C3 follow from ao and bo defined in i2. 17t by inter- 
changing v <-> A and f3 <-> a. In the Z' -plane, these two elliptic 
cylinders result in ellipses as in eq. d3.30t . with coefficients respec- 
tively given in eq. J3.3U and 



2 

C3 e\, 



dx 

d,2 = c 3 e 2 e3 + &3 (e? + e 2 ) 2 , 

d 3 = — 2c 3 eie 2 e 3 

d 4 = -2c 3 eie 2 (e? + e|) 2 (v z ,/g) 

d 5 = 2 e 3 (el + e| 



(3.40) 



1 ' 1 r c 3 e\ - b 3 (el + el)] (v z >/g), 



d 6 = (el + e 2 ) [(c 3 e 2 + 6363) (v z ,/gf - &3C3] . 

As before, A£' follows from the combination of the real zero 
points of the corresponding fourth order polynomials in eq. {333}, 
and of eq. Q32b in the case of maximum streaming around 
the long axis. For the outer long-axis tube orbits, S^ax = 
min[Sto P (A, fi, — /3), 5top(A, — a, v)]. The two elliptic cylinders 
are the one in eq. J3.29I ) and the second in eq. {335}, with the co- 
efficients of the corresponding ellipses in the Z'-plane are given in 
eq. {331} and eq. {336} (i = 2), respectively. 

The part of integral space accessible by box orbits is the region 
for which both v 2 , > at fi = — j3 and v 2 > at A = —a (Fig. [2}. 
Therefore, Sm&x = St op (A, jj,, v), and the two elliptic cylinders are 
the first in eq. {335} and the one in eq. {339}. The coefficients of 
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the corresponding ellipses in the Z'-plane are respectively those in 
eq. < f336t (i = 1) and in eq. < f3T40b . 



3.4 Gauss-Hermite moments 

We have seen that the line-of-sight velocity moments p,k(x' ,y') 
can be derived either via line-of-sight integration of the intrinsic ve- 
locity moments (eq. l3.10t or as moments of the LOSVD feq. 13 . lit . 
The lowest order line-of-sight velocity moments jtto, Mi and /12 pro- 
vide the surface mass density E, the mean line-of-sight velocity V 
and dispersion a by 

L = fio, V = — , and a — 5 , (3.41) 

Mo Mo 

all three as a function of (x' , y'). Whereas E, V and a can be mea- 
sured routinely, determinations of the higher order moments (/i3, 
Pi, ...) are more complicated. Spectroscopic observations of the 
integrated light of galaxies provide the LOSVD as function of po- 
sition on the sky plane. Unfortunately, the wings of the LOSVD 
become quickly dominated by the noise in the observations, and 
since the higher order moments significantly depend on the wings, 
their measurements can become very uncertain. Instead of these 
true higher-order moments, one often uses the Gauss-Hermite mo- 
ments (/13, hi, . . . ), which are much less sensitive to the wings of 
the LOSVD (van der Marel & Franx 1993; Gerhard 1993). 

There is no simple (analytic) relation between the true mo- 
ments and the Gauss-Hermite moments, including the lower order 
moments Egh, Vgh and <tgh (but see eq. 18 of van der Marel & 
Franx 1993 for approximate relations to lowest order in hz and hi). 
Nevertheless, we have shown that for Abel models the full LOSVD 
can be computed in a efficient way from eq. d3.27t , so that by fitting 
a Gauss-Hermite series to the resulting LOSVD, we can derive the 
Gauss-Hermite moments accurately, all as function of (x' , y'). 

Still, the calculation of the line-of-sight velocity moments 
through the intrinsic moments is useful, e.g., in case of investi- 
gating a range of viewing directions. The intrinsic moments have 
to be computed once, after which only a (numerical) integration 
along the line-of-sight is needed for each viewing direction. This 
is (much) faster than calculating the LOSVD separately at each di- 
rection. The higher order true moments can even be used to (nu- 
merically) determine the Gauss-Hermite moments. One way is to 
find the Gauss-Hermite LOSVD of which the true moments best- 
fit those from the Abel model. However, in practise this direct fit- 
ting of the true moments has several (numerical) problems. Because 
it is a non-linear minimisation problem, the convergence can take 
long and may result in a local instead of the global best-fit solution, 
possibly resulting in Gauss-Hermite moments that are significantly 
different from their true values. If, instead, we first (re)construct 
the LOSVD from the true moments by means of an Edgeworth 
expansion (see Appendix |C]( and then fit a Gauss-Hermite series, 
the Gauss-Hermite moments can be calculated accurately and effi- 
ciently. Evidently, once the viewing direction is known, it is more 
straightforward to compute the full LOSVD to derive the (higher- 
order) Gauss-Hermite moments. 

When we construct a galaxy model consisting of multiple 
Abel DF components (§ |2.4t , we cannot simply combine the cor- 
responding Gauss-Hermite moments in a linear way, because they 
are non-linear functions of the DF. Instead, we first add together the 



LOSVDs of the different DF componentjf], each multiplied with 
a constant weight, and then parameterise the resulting combined 
LOSVD as a Gauss-Hermite series. Because the mass included in 
each DF component is different, in order to obtain the mass frac- 
tions per DF component, we multiply the latter weights with the 
mass of the corresponding DF component divided by the total (lu- 
minous) mass. To change the sense of rotation of a rotating DF 
component (LR or SR), the corresponding observables do not have 
to be recomputed, as a change in the sign of the odd velocity mo- 
ments is sufficient. 



3.5 Surface brightness 

The surface brightness follows upon integration of the luminos- 
ity density along the line-of-sight. The luminosity density in turn 
is related to the mass density p* via the stellar mass-to-light ra- 
tio Mi,/L. With pi, the zeroth-order velocity moment of the DF 
(eq. l2.8t . the surface brightness follows as 

/oo 
(M,/L)"Vooo(x', y, z) dz' . (3.42) 
- 00 

In the special case when (Ah/L) does not change (e.g., due to 
variation in the underlying stellar populations) with position, we 
can take it out of the integral and SB = E/(A/*/L), where E is 
the surface mass density defined in eq. ( 13 .4 1 l l - 

In addition to the luminous matter, a galaxy may also contain 
dark matter. While in the outer parts of late-type galaxies the pres- 
ence of dark matter, as predicted by the cold dark matter paradigm 
for galaxy formation (e.g., Kauffmann & van den Bosch 2002), 
was demonstrated convincingly already more than two decades ago 
(e.g., van Albada et al. 1985), the proof in the outer parts of early- 
type galaxies remains uncertain, mainly due to a lack of kinematic 
constraints. As a consequence, in the outer parts of galaxies, com- 
monly a simple functional form for the dark matter distribution 
is assumed, often the universal profile from the CDM paradigm 
(Navarro, Frenk & White 1997). 

The dark matter distribution in the inner parts of galaxies is 
probably even more poorly understood (e.g., Primack 2004). For 
this reason, in current dynamical studies of the central parts of 
early-type galaxies, it is commonly assumed that both (M* /L) and 
the dark matter fraction are constant, i.e., mass follows light. In this 
case the surface brightness also follows from SB = Es/(M/L), 
where (M/L) is the (constant) total mass-to-light ratio and Eg the 
surface mass density, which after deprojection yields ps, the mass 
density related to the potential Vs via Poisson's equation $2A\ . In 
case of a Stackel potential i2.3\ . Es (and hence the surface bright- 
ness) has concentric isodensity contours that show no twist (e.g., 
Franx 1988). 



4 TRIAXIAL THREE-INTEGRAL GALAXY MODELS 

After choosing a Stackel potential, we investigate the shape of the 
density generated by the Abel DF components, and use these com- 
ponents to construct a triaxial galaxy model with three integrals of 
motion. 



Or, in case the LOSVD is not readily accessible, the true line-of-sight 
velocity moments, which are also linear functions of the DF. 
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4.1 Isochrone potential 

There are various choices for the potential that provide useful test 
models for comparison with the kinematics of triaxial elliptical 
galaxies (e.g., Arnold et al. 1994). One option is to consider the so- 
called perfect ellipsoid, for which Statler (1987) already computed 
numerical Schwarzschild models and Hunter & de Zeeuw (1992) 
investigated the maximum streaming thin orbit models. It has a 
density distribution stratified on similar concentric ellipsoids, but 
the potential function U (r) contains elliptic integrals, which slows 
down numerical calculations. An alternative is to consider the set 
of models introduced by de Zeeuw & Pfenniger (1988), which have 
nearly ellipsoidal density figures, and have a potential and density 
that are evaluated easily and swiftly. They are defined by the choice: 



U(t) = -GMVt(t + p), 

so that the triaxial Stackel potential has the elegant form 



V s {X,p,,v) = 



GM ( ^/Xp, + y/Jw + s/vX - p 



(Vx + y/p)(^/p + \/v)(sJv + VX) ' 



(4.1) 



(4.2) 



where we set GM = y/— 7 + \J—<x so that Vs(— a, — /?, —7) = 
— 1 in the centre. In the oblate axisymmetric limit this potential 
is that of the Kuzmin-Kutuzov (1962) models of Dejonghe & de 
Zeeuw (1988), and in the spherical limit it reduces to Henon's 
(1959) isochrone. For all these models, Vs = U[— a, — P,t] 
along the short z-axis is identical to the isochrone potential 
—GM/(t/t + \J—a). We therefore refer to models with U (r) of 
the form t4. It as isochrone models. Since the potential falls of as 
1/r at large radii, all these models have finite total mass. 

The expressions for the integrals of motion are given in i2.5l , 
where U[X,p.,u] — Vs and the third order divided difference 
U[\, fi, v, a] is given by the symmetric expressiorQ 



U[X,n, v, a] 



-GM-Vs(Vx + y/jI+Vt+Vt) 

(\fx + ^)(vfji+Vv)(Vv + Vv) ' 



(4.3) 



These triaxial isochrone models have the convenient property that 
the expressions for the potential and the integrals of motion contain 
only elementary functions of the (confocal ellipsoidal) coordinates 
and have no singularities. 

The same is true for the associated mass density ps, of which 
the expression is given in Appendix C of de Zeeuw & Pfenniger 
(1988), and a contour plot of ps in the (x, z) -plane is shown in 
their Fig. 2. These authors also derive the axis ratios of ps in the 
centre (their eq. C7) and at large radii (their eq. Cll), in terms of 
the axis ratios £ and £ of the confocal ellipsoidal coordinate system, 
defined as 



C = (-P)/(-a), e = (-7)/(-a). 



(4.4) 



Although ps becomes slightly rounder at larger radii, its axis ratios 
remain smaller than unity (for £ < £ < 1) because at large radii 
ps ~ 1 /r 4 in all directions. Characteristic values for the axis ratios 
can be obtained from the (normalised) moments of inertia along the 
principal axes of the density, 



/ x 2 p(x,0,0) dx 



(4.5) 



/ p(x, 0, 0) da; ' 
where the intermediate and short semi-axis length, b and c, of 



Substituting eq. A4.21 shows that U[X, p, v, a] is in fact fully symmetric: 

U[\,ti,i>,a] VA^iv4- x /^if(r + Vf<TA+ ^/a — /3( vA-f- y^T+ V^) 



the inertia ellipsoid follow from the long semi-axis length a by 
replacing x with y and z, and at the same time p(x, 0, 0) with 
p(0, y, 0) and p(0, 0, z), respectively. Taking for example £ = 0.8 
and £ = 0.64, the semi-axis lengths of the inertia ellipsoid result 
in the characteristic axis ratios bs/as = 0.88 and cs/as = 0.80 
for the density ps - The contours of the projected density are nearly 
elliptic with slowly varying axis ratios. 

For triaxial mass models with a Stackel potential Vs, de 
Zeeuw, Peletier & Franx (1986) have shown that the corresponding 
intrinsic mass density ps cannot fall off more rapidly than 1/r 4 , 
except along the short z-axis. All models in which ps falls off less 
rapidly than 1/r 4 become round at large radii. When ps ~ 1/r 4 , 
as is the case for, e.g., the above isochrone potential and the per- 
fect ellipsoid (e.g., de Zeeuw 1985a), the model remains triaxial at 
large radii. Moreover, mass models containing a linear combination 
of different Stackel potentials are possible as long as the associated 
confocal ellipsoidal coordinate systems share the same foci (e.g., de 
Zeeuw & Pfenniger 1988; Batsleer & Dejonghe 1994). This shows 
that, although we choose here a (single-component) isochrone po- 
tential, our method is capable of providing Abel models for a large 
range of Stackel potentials, with a similarly large range of shapes 
of the corresponding mass model. The same holds true for the lu- 
minous mass density, which we consider next. 



4.2 The shape of the luminous mass density 

Whereas the shape of the (total) mass density ps is fixed by the 
choice of the potential Vs, and £ and £ (eq. l4.4t . the shape of the 
(luminous) mass density p*, which is the zeroth order velocity mo- 
ment of the DF (eq. |2.8l l, also depends on the DF parameters w, u 
and S, and the type of component. For £ = 0.8 and £ = 0.64, we 
show in Fig. [4] for non-rotating DF components the characteristic 
(eq. 14. 5\ axis ratios of the corresponding density, as function of w 
and it. We have set 5 = 1, but the axis ratios depend only weakly 
on it, with p* becoming slightly flatter for increasing 5. The thick 
contours are drawn at the levels that correspond to the values of 
the characteristic axis ratios of ps, respectively bs/as = 0.88, 
cs/bs = 0.90 and cs/as = 0.80. These values are independent 
of w and u (as well as the other DF parameters). 

While the intermediate-over-long axis ratio b/a increases with 
increasing w, its value is only weakly dependent of it. By contrast, 
the short-over-intermediate axis ratio c/b mainly increases with in- 
creasing u. The short-over-long axis ratio c/a is the product of the 
previous two axis ratios and thus depends on both w and u. When 
both w and u are negative, the density p* has its long-axis along the 
x-axis and its short-axis along the z-axis, in the same way as the 
potential Vs and the associated density ps- Above certain positive 
values of either w or u, the axis ratios become larger than unity, 
which means that p* is no longer aligned with the underlying co- 
ordinate system in the same way as Vs and ps - For example, when 
(-a)tu = -0.5 and {-a)u = 0.5, b/a < 1 but c/b > 1, so that 
in this case p* has its short axis along the y-axis. 

A change in the sign of w and it has a strong effect on the ra- 
dial slope of p*. In Fig. [5] the radial profiles of p* along the princi- 
pal axes are shown for three combinations of w and it. The density 
is normalised to the central value po. The profiles along the y-axis 
(dotted curves) and along the z-axis (dashed curves) are arbitrarily 
offset vertically with respect to the profile along the a;-axis (solid 
curves) to enhance visualisation. The thin curves are the profiles of 
the (luminous) mass density p* for varying S, from S — (dark- 
est curve) to 8 — 4 (lightest curve), in unit steps. The thick black 
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Figure 4. The characteristic axis ratios b/a, c/b and c/a of the luminous mass density for a non-rotating Abel component, as function of the DF parameters 
w and u, while S = 1. The axis ratios of the confocal ellipsoidal coordinate system are £ = 0.8 and £ = 0.64, so that cf. j2.12t (— oi)w > —25/9 ~ —2.78 
and (— a)u < 625/144 pa 4.34. The thick contours are drawn at the levels that correspond to the characteristic axis ratios of the total mass density pg, 
associated with the underlying isochrone Stackel potential j4.2K respectively bg /ag = 0.88, cg/bg = 0.90 and cg/ag = 0.80. The intermediate-over-long 
axis ratio b/a depends mainly on w, the short-over-intermediate axis ratio c/b depends mainly on u, and c/a is the product of the previous two. 
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r/(-a) 1/2 r/(-a) 1/2 r/(-a) 1/2 

Figure 5. Principal axes profiles of the luminous mass density p + for a non-rotating Abel component, normalised to the central value p*,o- Each panel is 
for a different combination of the DF parameters w and u, while the grey scale indicates variation in S from zero (darkest curve) to four (lightest curve), in 
unity steps. The profiles along the j/-axis (dotted curves) and along the z-axis (dashed curves) are arbitrarily offset vertically with respect to the profile along 
the x-axis (solid curves) to enhance visualisation. The thick black curves show the profiles for the (total) mass density pg , associated with the underlying 
isochrone Stackel potential 14. 2t . with £ = 0.8 and § = 0.64. When the value of either w or u is positive (right panel), the profiles show a break at r ~ yj — a, 
so that these compact components may be used to represent kinematically decoupled components. 



curves show the profiles for the (total) mass density pg, which is 
independent of w, u and S. 

The profiles of p* steepen for increasing S and for increasing 
absolute values of w and u. In particular, when either w or u be- 
comes positive (right panel), the profiles suddenly become much 
steeper and drop to zero already at relatively small radii r ~ y/—a. 
The resulting Abel components are thus compact and, as we saw 
above, can be different in shape and orientation from the main body 
of the galaxy model. Therefore, they can be used to represent kine- 
matically decoupled components. When both w < and u < 
(left and middle panel), p* falls off much more gently and the Abel 
components cover a larger region. When w = u = (left panel), 



so that the DF only depends on energy, the profiles as well as the 
shape (Fig. [4) of p* can even be flatter than those of pg. However, 
already for small non-zero values of w and u, generally p* < ps 
everywhere in the galaxy model, and p* < ps in the outer parts. 
Although self-consistency p* = pg is only possible in the spher- 
ical case (for fixed values of w and u, see § \2.3l , one can choose 
the parameters w, u and 8 so that p* ~ ps- At the same time, hav- 
ing pi, < ps in the outer parts of the galaxy model, allows for a 
possible dark halo contribution. 

The shape of p* can furthermore change due to the additional 
contribution from long-axis rotating and short-axis rotating compo- 
nents. Although these components have no density along their ro- 
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tation axis, the behaviour of their overall shape as function of w, u 
and S is similar as for the corresponding non-rotating components. 

The above analysis shows that, given the triaxial isochrone 
potential ( 14.2b . we can use Abel components to construct a galaxy 
model with a realistic density. Depending on the choice of w, u 
and S, the galaxy model can contain compact (kinematic ally de- 
coupled) components and account for possible dark matter (in the 
outer parts). Furthermore, we show below that even with a small 
number DF components, enough kinematic variation is possible to 
mimic the two-dimensional kinematic maps of early-type galaxies 
provided by observations with current integral-field spectrographs. 
This means that we can construct simple but realistic galaxy models 
to test our Schwarzschild software (§[5land l6l4l l. 

4.3 A triaxial Abel model 

As before, we choose the isochrone Stackel potential ( 14.2b . we take 
£ = 0.8 and £ = 0.64 for the axis ratios of the coordinate system 
( |4.4t , resulting in a triaxiality parameter \2.6l of about T = 0.61, 
and we set the scale length \J—a — 10". Assuming a distance 
of D = 20 Mpc and a total mass of 10 11 Mq results in a central 
value for the potential Vo ~ 2.7 x 10 5 km 2 s -2 , which also sets 
the unit of velocity. We restrict the number of DF components to 
three, one of each type. For the first component of type NR we set 
W = u = — 0.5/(— q) and 8 = 1, so that the shape of the corre- 
sponding density is similar to that of ps, except in the outer parts 
where a steeper profile mimics the presence of dark matter (see 
Figs. E] and |5). For the second and third component, respectively 
of type LR and SR, we adopt the same parameters, expect that we 
take w = 0.5/(— a) and u — — 1.0/(— a) for the SR component, 
which therefore is more compact than the NR and LR component. 

We set the line-of-sight by choosing t? = 70° and ip = 30° 
for the viewing angles. After rotation over the misalignment an- 
gle tp = 101° eq. J3.1b , we compute for each DF component the 
LOSVD as a function of the positions on a rectangular grid on the 
sky plane, illustrated in Fig. [6] for five sky positions. By fitting a 
Gauss-Hermite series to each LOSVD, we obtain the maps of the 
mean line-of-sight surface mass density E, velocity V, dispersion 
a and higher-order Gauss-Hermite moments /13 and h&, shown in 
Fig. [7] The parameters of each DF component are given on the 
right. The NR component has zero (green) odd velocity moments. 
For the LR and SR component, the even velocity moments show 
a decrease in the centre, because these components have zero den- 
sity along respectively the intrinsic long and short axis. We add the 
LOSVDs of the NR, LR and SR components, weighted with mass 
fractions of respectively 80%, 12.5% and 7.5%, and fit a Gauss- 
Hermite series to obtain maps of E, V, a, /13 and hi. We convert 
E to the surface brightness by dividing by a constant stellar mass- 
to-light ratio of (A/*/L) = 4 M /L . 

To convert these 'perfect' kinematics to 'realistic' observa- 
tions, similar to those obtained with integral-field spectrographs 
such as SAURON (Bacon et al. 2001), we finally apply the following 
steps. We compute the kinematics on a rectangular grid consisting 
of 30 by 40 square pixels of l"in size. Using the adaptive spatial 
two-dimensional binning scheme of Cappellari & Copin (2003), we 
bin the pixels according to the criterion that each of the resulting 
(Voronoi) bins contains a minimum in signal-to-noise (S/N), which 
we take proportional to the square root of the surface brightness. 
For the mean errors in the kinematics we adopt the typical values 
of 7.5 kms -1 for V and a and 0.03 for /13 and hi in the kine- 
matics of a representative sample of early-type galaxies observed 
with SAURON (Emsellem et al. 2004). We then weigh these val- 



ues with the S/N in each bin to mimic the observed variation in 
measurement errors across the field. Finally, we use the computed 
measurement errors to (Gaussian) randomise the kinematic maps. 
In this way, we include the randomness that is always present in 
real observations. The resulting kinematic maps are shown in the 
top panels of Fig. [8] Because of the eight-fold symmetry of the 
triaxial model, the maps of the even (odd) velocity moments are 
always point-(anti)-symmetric, apart from the noise added. 



5 RECOVERY OF TRIAXIAL GALAXY MODELS 

We briefly describe our numerical implementation of Schwarz- 
schild's method in triaxial geometry (see vdB07 for a full descrip- 
tion), which we then use to fit the observables of the triaxial Abel 
model constructed in § 14.31 We investigate the recovery of the in- 
trinsic velocity moments and, through the distribution of the orbital 
mass weights, the recovery of the three-integral DF. 

5.1 Triaxial Schwarzschild models 

The first step is to infer the gravitational potential from the observed 
surface brightness. We do this by means of the Multi-Gaussian Ex- 
pansion method (MGE; e.g., Cappellari 2002), which allows for 
possible position angle twists and ellipticity variations in the sur- 
face brightness. For a given set of viewing angles {■&, ip, ip) (see 
§ 13. It . we deproject the surface brightness and we multiply it by a 
mass-to-light ratio (M/L) to get the intrinsic mass density, from 
which the gravitational potential then follows by solving Poisson's 
equation. We calculate orbits numerically in the resulting gravita- 
tional potential. 

To obtain a representative library of orbits, the integrals of mo- 
tion have to be sampled well. The energy can be sampled directly, 
but since the other integrals of motion are generally not known, 
we start, at a given energy, orbits from a polar grid in the (x, z)~ 
plane, which is crossed perpendicularly by all families of (regular) 
orbits. We restrict ourselves to the region in the first quadrant that 
is enclosed by the equipotential and the thin orbit curves to avoid 
duplication of the tube orbits. To have enough box orbits to support 
the triaxial shape, we also start orbits by dropping them from the 
equipotential surface (Schwarzschild 1979, 1993). 

Assigning a mass weight 7j to each orbit j from the library, 
we compute their combined properties and find the weighted su- 
perposition that best fits the observed surface brightness and (two- 
dimensional) kinematics. However, the resulting orbital weight dis- 
tribution may vary rapidly, and hence probably corresponds to an 
unrealistic DF. To obtain a smoothly varying DF, we both dither 
the orbits by considering a bundle of integrated orbits that were 
started close to each other, and we regularise when looking for the 
best-fit set of orbital weights by requiring them to vary smoothly 
between neighbouring orbits (in integral space). Finally, the best- 
fit Schwarzschild model follows from the minimum in the (Chi- 
squared) difference between (photometric and kinematic) observ- 
ables and the corresponding model predictions, weighted with the 
errors in the observables. 

5.2 Fit to observables of a triaxial Abel model 

In this case, the gravitational potential is known and given by the 
isochrone Stackel potential Vs eq. \4.ll . However, to closely sim- 
ulate the Schwarzschild modelling of real galaxies, we infer the 
potential from a deprojection of an MGE fit of the surface mass 
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Figure 6. Line-of-sight velocity distribution (LOSVD) at five different positions (x' , y') on sky-plane (in arcsec at the top of each column) of three different 
Abel components. The isochrone Stackel potential 14. 2t is used, with f = 0.8 and £ = 0.64 (T = 0.61), and scale length ^/—a = 10". The model is placed 
at a distance of D = 20 Mpc and the adopted viewing angles are ■& = 70° and tp = 30°. From top to bottom the LOSVDs of a non-rotating (NR), long-axis 
rotating (LR) and short-axis rotating (SR) Abel component are shown, with the corresponding DF parameters w, u and 6 given on the right. The height of 
each LOSVD is normalised to unity, and a (dashed) Gaussian distribution with zero mean and the same dispersion as the LOSVD is shown as a reference. 



density Es generated by Vs. The resulting potential reproduces 
Vs to high precision, with relative differences less than 10 -3 . We 
compute a library of orbits by sampling 21 energies E via a log- 
arithmic grid in radius from 1" to 123" that contains >99.9 per 
cent of the total mass. At each energy, we construct a uniform polar 
start space grid of 7 radii by 8 angles within the first quadrant of 
the (x, z)-plane and drop box orbits from a similar uniform polar 
grid on the equipotential surface in the first octant. This results in 
a total of 21x7x8x2 = 2352 starting positions, from each of 
which a bundle of 5 3 orbits are started. Taking into account the two 
senses of rotation of the tube orbits, this results in a total 441000 
orbits that are numerically integrated in the potential. 

We sum the velocities of each bundle of orbits in histograms 
with 401 bins, at a velocity resolution of 10 kms -1 . We fit the 
weighted sum of the velocity histograms to the intrinsic mass 
density p*, which we obtain from a deprojection of an MGE fit 
to the observed surface brightness, multiplied with the (constant) 
(M+/ L) — 4 Mq/Lq. Simultaneously, we fit the projected val- 
ues of the velocity histograms to the observed surface brightness 
and higher-order velocity moments. Finally, at the same time, we 
regularise the orbital weights in E and in the starting positions by 
minimising their second order derivatives. The strenght of the reg- 
ularisation is given by the a smoothening parameter (e.g., Cretton 
et al. 1999), which we set to A = 0.1 (see vdB07). 

From Fig. [8] it is clear that the (simulated) observables of the 
triaxial Abel model (top panels) are very well matched by the best- 
fit triaxial Schwarzschild model (bottom panels). The signature of 



the kinematically decoupled component in the maps of the mean 
line-of-sight velocity V and Gauss-Hermite moment /13 is accu- 
rately fitted, as well as the kinematics of the main body up to h,4 
within the added noise (§ 14. 3t . Below we investigate how well 
the intrinsic velocity moments as well as the three-integral DF — 
which are not (directly) fitted — are recovered. Here, we keep the 
mass-to-light ratio and the viewing angles fixed to the input val- 
ues of the triaxial Abel model (§ |4.3t , while in vdB07 we vary 
these global parameters to study how well Schwarzschild' s method 
is able to determine them. 

5.3 Intrinsic velocity moments 

We calculate the intrinsic first and second order velocity moments 
of the Schwarzschild model by combining the appropriate moments 
of the orbits that receive weight in the superposition, and investigate 
how well they compare with the intrinsic velocity moments of the 
Abel model. In general, there are three first (vt) and six second 
order velocity moments (v s vt) (s,t = x,y,z). Combining them 
yields the six dispersion components a st of the velocity dispersion 
tensor, where a% = (v s v t ) — (v s )(vt). 

We first consider the (x, z)-plane, as it is crossed perpendicu- 
larly by all four (major) orbit families. Because (v x ) = (v z ) — 
o~xy = o-yz — 0, we are left with (v y ) perpendicular to the 
(x, z) -plane as the only non- vanishing mean motion and a zx in 
the (x, «)-plane as the only non-vanishing cross-term. The average 
root-mean-square velocity dispersion ctrms is given by <r RMS = 
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Figure 7. Maps of the surface brightness (SB; in 10 3 Lq pc — 2 ), mean line-of-sight velocity V and dispersion a (both in km s" 1 ), and higher order Gauss- 
Hermite moments /13 and ?i4, of the same three Abel DF components as in Fig. [6] obtained by fitting a Gauss-Hermite series to the LOSVDs at each (pixel) 
position on the plane of the sky. The numerical artifacts at the edges of the compact SR component (third row) disappear when combined with components 
that extent over the full field-of-view (see e.g. the top row of Fig. [8). 



(<r xx + <Ty y + a zz )/3. The ratio (v y ) /<j-rms of ordered-over- 
random motion is a measure of the importance of rotation for the 
gravitational support of a galaxy. In Fig.[9j the colours represent the 
values of this ratio in the (a;, z)-plane, for the input triaxial Abel 
model (left panel) and for the best-fit triaxial Schwarzschild model 
(right panel). 

In a Stackel potential the axes of the velocity ellipsoid are 
aligned with the confocal ellipsoidal coordinate system (e.g., Ed- 
dington 1915; van de Ven et al. 2003). As a result, one of the axes 
of the velocity ellipsoid is perpendicular to the (x, z) -plane, with 
semi-axis length a yv . The other two axes lie in the (x, z)-plane and 
have semi-axis lengths given by 



2 . 2 
+ cr 2z 



\2 1 4 
+ <7xz 



(5.1) 



The ellipses overplotted in Fig.|9]show the corresponding cross sec- 
tions of the velocity ellipsoid with the (x, z)-plane. The flattening 
of the ellipses is thus given by the ratio a-/a+, while the angle 
8 XZ of the major-axis with respect to the a>axis is given b}Q 



tan(26U) = 2a 2 xz /(a : 



.)• 



(5.2) 



In addition, the cross on top of each ellipse represents the ratio 



a yy /a+, i.e., the (relative) size of the velocity ellipsoid in the per- 
pendicular direction. For an isotropic velocity distribution the el- 
lipses become circles and the crosses fill the circles. Finally, the 
black curves are contours of constant mass density in steps of one 
magnitude. 

The density of the triaxial Abel model (solid curve) is well 
fitted by the triaxial Schwarzschild model (dashed curve), with a 
(biweighffj) mean fractional difference below 1 %. In both the Abel 
model and the fitted Schwarzschild model the value of (v y ) / ctrms 
is relatively low, with a mean value ~ 0.14, indicating that gravi- 
tational support is mainly due to random motion. Still, the average 
rotation of the long-axis tube orbits (with (v y ) < 0) due to the 
maximum streaming LR component in the input Abel model, as 
well as, the opposite maximum streaming of the (compact) short- 
axis rotating component are clearly visible, and well recovered by 
the best-fit Schwarzschild model. The average absolute difference 
in both (v y ) and ctrms is below 6 km s -1 , and thus well within the 
typical error of 7.5 km s _1 assigned to the simulated mean line-of- 
sight velocity V and velocity dispersion a of the Abel model (see 
W.3t . The corresponding uncertainty in (i^/orms is ~ 0.03. 

We see in Fig.|9]that, at larger radii, the ellipses become more 
radially elongated and the relative size of the crosses decreases in 
the radial direction, but they stil fill the ellipses in the angular direc- 



In case of alignment with the confocal ellipsoidal coordinate system, 
this angle is given by the tangent to the curves of constant u), i.e., 
tan X z = (z/x)(X + ct)/(A + 7), which indicates approaching align- 
ment with the polar coordinate system at large radius. 



8 The biweight mean (e.g., Andrews et al. 1972; Beers, Flynn & Gebhardt 
1990) is robust estimators for a broad range of non-Gaussian underlying 
populations and is less sensitive to outliers than other moment estimators. 
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Figure 8. Kinematic maps for a triaxial Abel model (top row) and converted to observables with realistic measurement errors added (middle row; see §[43}, 
and for the best-fit triaxial Schwarzschild model (bottom row; see § 15. It . From left to right: mean line-of-sight velocity V and dispersion a (both in km s — 1 ), 
and Gauss-Hermite moments hi and h,4. Isophotes of the surface brightness of the Abel model are overplotted in each map. At the right side of each map, the 
(linear) scale of the corresponding kinematics is indicated by the colour bar, and the limits are given below. 



tion. This implies a velocity distribution that becomes increasingly 
radially anisotropic outwards, but remains close to isotropic in the 
tangenetial direction everywhere. This shape and orientation of the 
velocity ellipsoid in the input Abel model is well reproduced by the 
best-fit Schwarzschild model, with only a (mild) underestimation 
of the radial anisotropy towards the z-axis. This is likely the re- 
sult of numerical difficulties due to the small number of (sampled) 
long-axis tube orbits that contribute in this region. The absolute 
difference in the semi-axis lengths a+, <j_ and a yy of the velocity 
ellipsoid is on average ~ 8kms _1 . This uncertainty includes both 
deviations in shape and orientation of the velocity ellipsoid, and is 
wihtin the expected range due to the errors in the simulated kine- 
matics. The corresponding axis ratios a-/a+ and cr yy /<j+ of the 
velocity ellipsoid are on average recovered within ~ 5 %. 



Away from the (a;, z)-plane, the average fractional differ- 
ence in the density between the input Abel model and the best-fit 
Schwarzschild model stays below 1%. Fig. 1101 compares the in- 
trinsic moments of the input triaxial Abel model (top) with those 
of the best-fit triaxial Schwarzschild model (bottom) in three di- 
mensions. The first column shows the (amplitude) of the streaming 
motion u str , given by v% T = v \ + v y + v%, and normalised by 
o"rms- These quantities are computed on a polar grid (r,9,(f>) in 
the first octant. The (logarithmic) sampling of the radius r is indi- 
cated by the black dots between the top and bottom panels, while 
each row is for a different polar angle 6 as indicated on the right, 
and the colours represent the (linear) change in azimuthal angle (f>. 
The limit (f> — 0° (black curves) corresponds to the (:r, z)-plane 
discussed above. The resulting ordered-over-random motion V/cr 
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Figure 9. The colours represent the mean motion (vy) perpendicular to the (x, z)-plane, normalised by ctrms (excluding the axes to avoid numerical 
problems), for the input triaxial Abel model (left) and for the best-fit triaxial Schwarzschild model (right). The ellipses are cross sections of the velocity 
ellipsoid with the (a;, z)-plane and the crosses represent the (relative) size of the velocity ellipsoid in the perpendicular (j/-axis) direction. The black curves 
are contours of constant mass density in steps of one magnitude, for the input Abel model (solid) and for the fitted Schwarzschild model (dashed). See § 15.31 
for details. 



is well recovered by the Schwarzschild model, apart from the upper 
panel, which is likely the result of the above mentioned numerical 
difficulties close to the z-axis. Overall, the average absolute differ- 
ence in both v BtI and urms is below 6 km s _1 and the uncertainty 

in v s tr/crRMS is ~ 0.03. 



The second and third column of Fig.[l0]show respectively the 
intermediate-over-major at,/ a a and minor-over-major a c /a a axis 
ratios of the velocity ellipsoid. The velocity ellipsoid of the triax- 
ial Abel model is aligned with the confocal ellipsoidal coordinate 
system, so that its semi-axis lenghts a a > ut > a c follow directly 
from a 2 — {v 2 } — {vt}' 2 with r = A, /i, v. In general, this is not 
the case for the triaxial Schwarzschild model, and instead we diago- 
nalize the (symmetric) velocity dispersion tensor with components 
{a st) (s, t — x,y, z). As before, the axis ratios of the velocity ellip- 
soid are quite well recovered by the best-fit Schwarzschild model, 
except towards the z-axis (upper panels) where it underestimates 
the anisotropy in the velocity distribution of the input Abel model. 
Similarly, away from the (a;, z)-plane {<j> = 0°, black curves), the 
Schwarzschild model increasingly overestimates the a^ja a ratio, 
while the a c /a a remains well reproduced. It is plausible that the 
recovery in the (x, z)-plane is better, because it is optimally sam- 
pled as starting space for the numerical orbit calculations, and it 
is crossed perpendicularly by all four major orbit families. Never- 
theless, the absolute difference in a a , <Jb and a c between the in- 
put Abel model and the best-fit Schwarzschild model is on average 
~ 9kms _1 . The axis ratios at/a a and a c /a a are on average re- 
covered within ~ 6 %. 



5.4 Three-integral distribution function 

The fitted triaxial Schwarzschild model results in a mass weight 7 
per orbit. These mass weights are a function of the three integrals 
of motion (E, l2,Is)- In general, only the energy is exact, but for 
a separable potential I2 and J3 are also known explicitly and given 
by eq. d2.5t . The orbital mass weights follow from the DF by in- 
tegrating f(E,l2,H) over the part of phase-space (x,v) that is 
accesible by the orbit. Since each orbit is a (unique) delta-function 
in integral-space, the resulting orbital mass weights are in principle 
zero. However, as described in § 15.11 and § 15.21 final orbits con- 
sists each of a bundle of 125 orbits started closely to each other and 
their assigned mass weights are required to vary smoothly between 
neighbouring orbits. 

To estimate the orbital mass weights from the input triaxial 
Abel model, we divide the integral-space in finite cells and link 
each cell to the orbit that corresponds to its centroid. The corre- 
sponding mass weights then follow from 

j(E,l2,I 3 )=J f jf{E,I 2 ,h) AV(E,I 2 ,h) dEdhdh, (5.3) 



where 



AV(E,I 2 ,I 3 ) = JJJ 



d(v x 



d{E,h,h) 



dx dy dz, 



(5.4) 



with SI the volume in configuration space accessible by the orbit. 
The multi-component DF of the input triaxial Abel model consists 
of basis functions defined in eq. J2.24b , with the DF parameters and 
weights per component given in § 14.31 Below, we first calculate A V 
and the cell in integral space, and then return to the comparison of 
the orbital mass weights. 
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Figure 10. Intrinsic velocity moments in three dimensions for the input triaxial Abel model (top) and for the best-fit triaxial Schwarzschild model (bottom). 
The first column shows the (amplitude) of the streaming motion i%tr , normalised by ctrms ■ The second and third column show the axis ratios of the velocity 
ellipsoid, where a a , cr b and cr c are respectively the semi-lengths of the major, intermediate and minor principal axes. These quantities are computed on a 
polar grid (r, 8, <j>) in the first octant. The (logarithmic) sampling of the radius r (in arcsec) is indicated by the black dots between the top and bottom panels. 
Each row is for a different polar angle 9 (in degrees) as indicated on the right, with the top panel close to the z-axis and the bottom panel close to the (x, y)- 
plane. The colours represent the (linear) change in azimuthal angle (p (in degrees), with limits 0° and 90° corresponding to the (x, z)-plane and (y, z)-plane, 
respectively. See ij |5.3l for details. 
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Figure 11. Three quantities involved in the calculation of the orbital mass weights for a triaxial Abel model with an isochrone potential. For a given energy 
E, in each panel, the values of the second and third integral of motion, I2 and 73, indicated by the symbols, correspond to the orbital starting position and 
velocities in the triaxial Schwarzschild model that is fitted to the observables of this triaxial Abel model. The solid curves, calculated with the expressions 
in § 15.4.21 bound and separate the regions of the box (B) orbits, inner (I) and outer (O) long-axis tube orbits and short-axis (S) tube orbits. The circles refer 
to orbits started in the (x, z)-plane and the triangles represent the additional set of orbits dropped from the equipotential surface (see § [3TTJ. The latter box 
orbits may overlap with those started from the (x, z)-plane. The colours in the left panel indicate the value of the DF f(E, I2, 13) for each orbit in units 
of Mq (km/s) -3 pc -3 , with the (linear) scale given by the vertical bar on the right. The colours in the middle panel represent the values of AV(E, I2, 13) 
defined in eq. (5A\ in units of (km/s)- 3 pc -1 . The area of each Voronoi bin in the right panel, multiplied by the range in energy E, approximates the cell 
AEA(l2, 13) in integral space for each orbit. The product of these three quantities yields an estimate of the mass weight "f(E, l2,Iz) for each orbit. 



5.4.1 Integral over configuration-space 

The expression for AV(E, 12,^3) of a single orbit in a triaxial 
Stackel potential can be deduced from the relations in § 7.1 of de 
Zeeuw (1985a). It is given by 



AV(E,I 2 ,h) = (7- 



{\ - - v){v - A) 
ffl(A) a(/x) a(v) 



X 



E-V ctt (\)] [E-Vesifj,)] [E-V eS (v)] 



where a(r), r = A, fj,, is defined as 

a(r) = (t + c*)(t + /3)(t + 7), 



dXdfidv, (5.5) 



(5.6) 



the effective potential V c s as 

h h 



Kff(T) 



+ 



+ U[t, -a, -7], 



(5.7) 



r + a t + 7 

and Q is the volume in configuration space accessible by the orbit 
in the triaxial separable potential that obeys (E, l2,Iz). The last 
term in eq. \5.1\ is equal to the Stackel potential l |2.3t along the 
intermediate j/-axis. 

Because of the separability of the equations of motion, each 
orbit in a triaxial separable potential can be considered as a sum 
of three independent motions. Each of these one-dimensional mo- 
tions is either an oscillation or rotation in one of the three confocal 
ellipsoidal coordinates (A, fi, v), such that the configuration space 
volume Q is bounded by the corresponding coordinate surfaces. 
The values of (A, fi, v) that correspond to these bounding surfaces 
can be found from Table [Tj for the four families of regular orbits: 
boxes (B), inner (I) and outer (O) long-axis tubes, and short-axis (S) 
tubes. Whereas a, f3 and 7 are the limits on (A, n, v) set by the foci 
of the confocal ellipsoidal coordinate system, the other limits are 
the solutions of E = V c b(t) (see Fig. 7 of de Zeeuw 1985a). In the 



case of the triaxial isochrone Stackel potential l |4.2t , we can write 
this equation as a fourth-order polynomial in \pr. The solutions are 
then the squares of three of the four roots of this polynomial (the 
fourth root is always negative). 

For each orbit in our Schwarzschild model, we compute 
(E, l2,H) by substituting the starting position and velocities of the 
orbit into the expressions {23} • From the value of E and the sign 
of I2 (while always 73 > 0), we determine to which orbit family it 
belongs. The corresponding configuration space volume Q. is then 
given by the boundaries for A, fj, and v in the last three columns of 
TableQ] The value of AV(E, I2, 13) follows by numerical evalua- 
tion of the right-hand side of eq. \5.5l . 

The integrand in eq. l |5.5l > contains singularities at the integra- 
tion limits, which can be easily removed for a triaxial isochrone 
potential. We write the integrand completely in terms of (-y/cr ± 
y/r) 1 / 2 , where a, r — A, fi, v or a constant value. Suppose now 
that the integral over A ranges from Ao to Ai and the terms (\[\ — 
^Ao") 1/2 and(VAT-v / A) 1/2 appear in the denominator. The sub- 
stitution \f\ — \/Ao + (\/A7 — %/Ao) sin 2 V then removes both 
singularities since 



5.4.2 Cell in integral space 

We approximate the triple integration over the cell in integral space 
in eq. l |5~3l l by the volume AEA(l2,h). Here AE is the (log- 
arithmic) range in E between subsequent sets of orbits at differ- 
ent energies (see § I5.lt . with limits given by the central poten- 
tial and E = 0. Because we do not directly sample I2 and ^3 
in our implementation of Schwarzschild's method, as their expres- 
sions are in general unknown, we cannot directly calculate the area 
A (I2, h). Instead, we compute the Voronoi diagram of the points 
in the (I2, /3)-plane that correspond to the starting position and ve- 
locities of each orbit, at a given energy E. An example is given in 
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Table 1. Configuration space volume Q accessible by the four families of regular orbits. 



orbit I2 E X fx v 



Box orbits 


< 


V eff (-/3) . 


. 




— a . 


• Vax 


-p . 


■ /^max 


-7 ■ 


■ Z-'max 


Inner long-axis tube orbits 


< 


min [y efI (/i)] . 


• V off (- 


-P) 


—a . 


- Amax 


Mmin ■ 


- ^max 


-7 • 


• -P 


Outer long-axis tube orbits 


> 


min[V eff (A)] . 


• V cB (- 


-0) 


■^min ■ 


• Vax 


A*min ■ 


. — OL 


-7 • 


■ -P 


Short-axis tube orbits 


> 


max{V cff (-/3),min[V eff (A)]} . 


. 




■^min ■ 


• Amax 


-P ■ 


. —OL 


-7 ■ 


■ ^max 



the right panel of Fig. QT| The area of the Voronoi bins approxi- 
mates the area A(/2, ^3) for each orbit. 

The four families of regular orbits are separated by two lines 
that follow from I2 = and E = Kff ' (— /3). The latter provides the 
part of the boundary on I2 and ^3 for the box orbits. The remainder 
of this boundary is given by the positivity constraint on J3 and by 
the solution of (cf. eqs 64 and 65 of de Zeeuw 1985a) 



E = Vcb(ko) and 



dVeff(«0 



(Ik 



0, 



ko > -p. 



(5.8) 



Substituting V e g from eq. J5.7t and using dU[r, —a, — j]/dr = 
U[t, t, —a, —7], we find the solution 



(kq + ay 
(a -7) 



{e-u[- 



a, k , K \ 



(5.9) 



and similarly for 73 by interchanging a <-> 7. For —p < kq < — a, 
the solution describes the boundary curve for which I2 < and 
corresponds to the thin I tube orbits. For kq > —a, we find the 
boundary curve for which I2 > 0, corresponding to the thin O and 
S tube orbits. 

There are limits on the values of Ko depending on the value 
of E, and sometimes there are no valid solutions for no, which 
implies that only box orbits contribute at that energy. These limits 
can be obtained from the thin orbit curves in the (a?, z)-plane. With 
y = v x = v z = 0, the expressions d2.5b for the integrals of motion 
reduce in this plane to 



E 
h 
h 



^Vy + U[\, K 



Pi 

{±v 2 y + (a-p)U[\,K,-p,-a\}, 
{± v * + ^ - P)U[\, K ,-p,-~,]}, 



(5.10) 



with —7 < k < —a replacing [i and v respectively above and 
below the focal curve given by z 2 /(■y — /?) — x 2 /(P — a) = 1. 
Next, we substitute the expression for E in those for I2 and ^3 and 
we use that (r + P)U[X, k, t, -0] = U[X, k, r] - U[X, k, —0\, 
respectively for r = — a and r = —7. We find that the thin orbit 
curves follow by solving I2 — and thus E — U[X, ko, ko] for 
I tubes, and I3 = and thus E — U[ko, ko, k], with k = fj, for 
O tubes and k — v for S tubes. In general these equations have 
to be solved numerically, but in the case of the triaxial isochrone 
potential l |4.2t , they reduce to a second order polynomial in y^^o 
and the solutions simply follow from the roots of the polynomial. 



5.4.3 Orbital mass weight distribution 

Once we have computed for each orbit the DF f(E,l2,l3), 
AV(E,I 2 ,Is) and the cell AEA(I 2 ,h) in integral space 
(Fig. II It . its (approximate) mass weight y(E, I2, 13) follows by 
multiplication of these three quantities. As before, the choice of 
maximum streaming for the (LR and SR) rotating components re- 



duces the accessible integral space, and thus also the corresponding 
orbital mass weights, by a factor two. 

The resulting orbital mass weight distribution of the input tri- 
axial Abel model is shown in the top panels of Fig. [T2] and that 
of the fitted triaxial Schwarzschild model in the bottom panels. The 
energy E increases from left to right, which corresponds to increas- 
ing distance from the centre as is indicated by the radius Re (in 
arcsec) at the top of each panel. For this representative radius we 
use the radius of the corresponding thin (S) tube orbit on the a;-axis. 
The values of I2 and ^3, on the horizontal and vertical axes respec- 
tively, are both normalised with respect to their maximum ampli- 
tude at the given energy. In each panel, the mass weight values are 
normalised with respect to the maximum in that panel. Between the 
two rows of panels, the fraction of the summed values in each panel 
with respect to the total mass weight in all panels is given (in %). 

The panels with Re < 40" are best constrained by the kine- 
matic observables. This takes into account that even orbits that ex- 
tend beyond the maximum radius covered by the observables can 
contribute significantly at smaller radii. In these panels, the main 
features of the orbital mass weight distribution of the triaxial Abel 
model are recovered. In the outer parts the Schwarzschild model 
is still constrained by the mass model, which extends to a radius 
of about 100", but the orbital mass weight distribution deviates 
from that of the input Abel model due to the lack of kinematic con- 
straints. A point-by-point comparison yields an average fractional 
error of ~ 50 %, and if we consider in each panel the mass weights 
above the mean value, which together contribute more than half 
of the total mass, the fractional error decreases to ~ 30 %. How- 
ever, this way of quantifying the recovery is (somewhat) mislead- 
ing since the relatively large fractional errors are at least partially 
caused by the strong peaks in the orbital mass weight distribution. 
For example, if in the input Abel model a certain orbit gets a sig- 
nificant weight, but in the Schwarzschild model, due to numerical 
uncertainties, this weight is assigned to a neighboring orbit with 
a (slightly) different value of ^3, the relative error at each of the 
corresponding points in the integral space can be very large. 

Henceforth, we show in Fig. [O] the orbital mass weights as 
function of each of the three integrals of motion separately by col- 
lapsing the cube in (E, 12,13) in the remaining two dimensions. 
We again use Re as a representative radius for E (first panel), but 
since the (range of) values for I2 and ^3 change with E (see also 
Fig.ll2t. we use their index in the cube instead. In addition to the 
total distribution, we also show the contribution of the three differ- 
ent NR, LR and SR components separately, as well as for the latter 
two rotating components the contributions from the two directions 
of rotation by making the mass weights for one of the directions 
negative. Since the input triaxial Abel model (diamonds connected 
by solid curves) is constructed with maximum streaming in one of 
the two directions for both the LR and SR component (see £j 14. 3) , 
the opposite direction in both cases has zero mass weight. This 
is nicely reproduced by the best-fit triaxial Schwarzschild model 
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Figure 12. The orbital mass weight distribution for the input triaxial Abel model (top) and for the best-fit triaxial Schwarzschild model (bottom). From left 
to right the energy increases, corresponding to increasing distance from the centre, indicated by the radius Re (in arcsec) of the thin short-axis tube orbit on 
the x-axis. The vertical and horizontal axes represent respectively the second and third integral of motion, I2 and 73, normalised by their maximum amplitude 
(for given E). In each panel, the colours represent the mass weights, normalised with respect to the maximum in that panel, and with the (linear) scale given 
by the vertical bars on the right. Between the two rows of panels, the fraction (in %) of the included mass with respect to the total mass is indicated. 



(crosses connected by dotted curves) in which ~ 2 % of the total 
mass, or ~ 10 % of the mass of the LR and SR components, is 
wrongly assigned to the opposite direction. Keeping in mind that 
the orbital mass weights itself are not directly fitted and that the 
typical velocity error of 7.5 km s _1 is more than 10 % of the max- 
imum in the simulated velocity field (see ^ 14.3b . these per centages 
are (well) within the expected uncertainties. 

From the first panel of Fig. [T3] we see that mass as function 
of E is well recovered, even in the outer parts where (nearly) all 
the constraints come only from the mass model. The average abso- 
lute difference is ~ 0.7 %. Whereas for E the constraints provided 
by the mass model already seem sufficient, for I2 and J3 the kine- 
matic constraints are essential. Not suprisingly, we then also see 
that the recovery is less good with an average absolute difference 
of ~ 1.9 % in I2 and ~ 1.0 % in I2. The main contribution is from 
the NR component, while the two rotating components seem to bet- 
ter recovered. 



6 AXISYMMETRIC THREE-INTEGRAL GALAXY 
MODELS 

We now consider three-integral galaxy models in the axisymmet- 
ric limit. As we have seen in the Introduction (Section QJ, various 
groups have successfully developed independent axisymmetric im- 
plementations of Schwarzschild's method and verified their codes 
in a number of ways. The published tests to recover a known (ana- 
lytical) input model have been limited to spherical geometry or to 
an axisymmetric DF that is a function of the two integrals of motion 
E and L z only. 

Here, we present the velocity moments of the three-integral 



Abel DF in the axisymmetric limit and we choose again the 
isochrone form in eq. ( 14.11 ) for the Stackel potential. The proper- 
ties of the resulting three-integral Kuzmin-Kutuzov models can be 
expressed explicitly in cylindrical coordinates. We construct an ax- 
isymmetric oblate Abel model and fit Schwarzschild models to the 
resulting observables to test how well the axisymmetric implemen- 
tation of Schwarzschild's method, as presented in Cappellari et al. 
(2006), recovers the intrinsic velocity moments as well as the three- 
integral DF. 



6.1 Velocity moments and line-of-sight velocity distribution 

When two of the three constants a, (3 or 7 are equal, the confocal 
ellipsoidal coordinates (A, /1, v) reduce to spheroidal coordinates 
and the triaxial Stackel potential l |2,3t becomes axisymmetric. 



6.1.1 Oblate axisymmetric model 

When j3 — a 7^ 7 (triaxiality parameter T = 0), we cannot use \x 
as a coordinate and replace it by the azimuthal angle cj>, defined as 
tan0 = y/x. The relation between (A, <f>, v) and the usual cylin- 
drical coordinates (R, <f>, z) is given by 



(X + a)(u + a) 



(A + 7)^ + 7) 



(6.1) 



The Stackel potential Vs(A, v) — U[X, —a, v] is oblate axisym- 
metric. The corresponding integrals of motion follow by substitu- 
tion of /1 = —/3 = —a in the expressions J2.5b . so that the second 
integral of motion reduces to I2 = \L%. 

With the choice J2.9I ) for the DF, the expression for ve- 
locity moments fJ,imn(X,v) is that of the triaxial case given in 
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Figure 13. The orbital mass weights (in % of the total mass) for the input triaxial Abel model (diamonds connected by solid curves) and for the best-fit triaxial 
Schwarzschild model (crosses connected by dotted curves), as function of each of the three integrals of motion. These 'projections' of the three-dimensional 
orbital mass weight distribution shown in Fig. ll2l are obtained by collapsing the cube in (E, l2,I:j) in two dimensions. As before, we represent the energy E 
in the first panel by the radius Re (in arcsec) of the thin short-axis tube orbit on the x-axis. For the second and third integral of motion, I2 and 73, we use 
the index in the cube, since the (range of) their values changes with E. The total distribution (black colour) is split into contributions from the non-rotating 
(NR; red), long-axis rotating (LR; green) and short-axis rotating (SR; blue) components. Moreover, for each rotating component the contributions from the 
two directions of rotation are separated by making the mass weights for one of the directions negative. 



eq. ( 12.10b . but with /i = — j3 = —a. From Fig. [T] we see that 
the lower limit on w vanishes. For the NR type of components, 
Smax = Stop (A, fi, —7) and the corresponding velocity moments 
/U(JrO,(A, v) vanish when either I, m or n is odd. Because the only 
family of orbits that exists are the short-axis tube orbits, we can 
introduce net rotation (around the short 2-axis) by setting the DF 
to zero for L z < 0, so that fj,f^ n (X, v) = §/4^(A, v). These SR 
velocity moments vanish when either I or n is odd, but are non-zero 
if m is odd. They should be multiplied with (— l) m for maximum 
streaming in the opposite direction. By choosing different weights 
for both senses of rotation, we can control the direction and the 
amount of streaming motion. 

In the conversion to observables described in § [3] the ma- 
trix Q, which transforms the velocity components (v\,v^,Vv) to 
(v x ,v y ,v z ), reduces to 



Q 



'A cos < 
A sin ( 
v B 



— sm ( 
cos <f) 





where A and B are defined as 

2 _ (A + 7 )0 + a) 



.4 



B 2 = (A + a)(i/ + 7 ) 



(6.2) 



(6.3) 



(A-v)(a-7)' " (A-i/)(7-a)' 

Because of the symmetry around the short-axis, the azimuthal 
viewing angle ip looses its meaning and the misalignment angle 
-0 = 0°. We are left with only the polar viewing angle which is 
commonly referred to as the inclination i, with i = 0° face-on and 
i = 90° edge-on viewing. As a consequence, the projection matrix 
P is a function of i only and follows by substituting 7? = i and 
ip — in eq. J3,4t . The rotation matrix R in eq. i3.5\ reduces to 
the identity matrix, so that M = PSQ. 

The expression for the LOSVD follows from that of the tri- 
axial case in eq. J3.27t by substituting /1 = — f3 — —a. For the 
NR components, again A£^r = 2-7T and the simplified expression 
< !3.28t holds in case of a DF basis function as defined in eq. J2.24| >. 
To introduce net rotation, we require that (w M =) v$ > as in 



§ 13.3.41 which yields SR components with maximum streaming. As 
illustrated in the right panel of Fig. [3] A£sr is the length of the part 
of the circle between the intersections £± = 2 arctan(ti± ) with the 
line (with u± given in eq. l3.38b . and which is on the correct side of 
the line in eq. ( 13.37b . This is again similar to SR components in the 
triaxial case, but without the restriction to stay within the ellipses. 



6.1.2 Prolate axisymmetric model 

When P = 7 7^ a (T = 1), we replace the coordinate v by 
the angle x, defined as tanx = z/y. The resulting coordinates 
(A, fi, x) follow from the above coordinates (A, (j>, v) by taking 
v — > fi, tj> — > x» an d 7 — > a — *• p. The Stackel potential 
Vs(X, n) = U[X, n, —7] is now prolate axisymmetric. By substi- 
tuting v — —j3 — —7 in eqs < |2.5 b and J2.10t , we obtain the ex- 
pressions respectively for the integrals of motion (with J3 = \B 2 X ) 
and for the intrinsic velocity moments /i; m „(A, /i). From Fig. Q] we 
see that now the upper limit on u vanishes. For the NR components, 
5 max = Stop (A, fi, —7), and since we only have the long-axis tube 
orbits, we can introduce net rotation (around the a;-axis) by setting 
the DF to zero for L x < 0, so that fi^ n (X, fi) — |/i^ n (A, /x). 
These LR velocity moments vanish if either I or m is odd and mul- 
tiplication with (—1)" yields net rotation in the opposite direction. 

The matrix Q, which transforms (vx, v^, v x ) to (v x ,v y , v z ), 
in this case reduces to 



C -D \ 

Q = ID cos x C cos x — sin x 
\Dsinx Csinx cosxy 



where C and D are given by 



(A + /3)( m + q) 
(A-/i)(a-/3)' 



(A + <»)(/* + /?) 
(A -/*)(/?- a)" 



(6.4) 



(6.5) 



In the projection matrix P in eq. l |3.4| l. we substitute # = 7r/2 — i 
and <p = 0, so that for inclination 2 = 0° and i — 90°, we are 
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respectively viewing the prolate mass model end-on and side-on. 
In the rotation matrix R we take ip = 90° to align the projected 
major axis horizontally. The expression for the LOSVD follows 
from eq. d3.27t by substituting v = — /3 = —7, and by requir- 
ing {v v =) v x — we obtain LR components with maximum 
streaming. As for the oblate case and illustrated in the left panel 
of Fig. [3] A£g R is the length of the circle part between the angles 
£± = 2arctan(u±) (with u± given in eq. I3.34t which is on the 
correct side of the line in eq. ( 13 .32b . 



6.2 Kuzmin-Kutuzov potential 

In the axisymmetric limit, the form d4.lt for U(t) results in the 
Kuzmin-Kutuzov (1962) potential. We give the properties relevant 
for our analysis, while further details can be found in Dejonghe & 
de Zeeuw (1988), including expressions and plots of the mass den- 
sity ps, its axis ratios, and the two-integral DF f(E, L 2 ) consistent 
with ps [see also Batsleer & Dejonghe (1993), who also corrected 
a typographical error in f(E, L 2 )]. 

When j3 = a, the oblate axisymmetric potential Vs(A, v) = 
U[X,—a,u] and the third order divided difference U[X, — a, v, a], 
which both appear in the expressions for the integral of motions 
d2.5t . have the simple forms 

-GM 



V s {X,u) = 



U[X, —a, v, a] 



VX + y/v' 



GM 



(VA + y/v)(VX + xfo)(yfv + y/a) 



(6.6) 



(6.7) 



where again GM = ^/— 7 + V~ a > so mat Vs = — 1 in the centre. 
By means of the relations 

A + v — R 2 + z 2 — a — 7, Xv — ccy — -yR 2 — az 2 , (6.8) 

and (y/X + y/v) 2 = A + i/ + 2^Vand (a/A + y/o)(y/u + y/a) = 
\f~Xv + y/a(y/X + s/v) + a, we can write the potential and in- 
tegrals of motion explicitly as elementary expressions in the usual 
cylindrical coordinates. 

When j3 = 7, the prolate potential Vs(X,fi) = U[X,fi, —7] 
and the third order divided difference U[X, fi, —7, a] follow respec- 
tively from d6,6t and d6,7t by replacing v by /i. 



6.3 An axisymmetric Abel model 

The above constructed triaxial Abel model (§ 14. 3t transforms into 
an oblate axisymmetric Abel model if we let £ approach unity, 
while keeping £ = 0.64 fixed. Similar to the triaxial case, the DF 
contains a NR component with the same parameters, u = w = 
— 0.5/(— a) and 5 = 1, but we exclude the LR component since 
long-axis tube orbits do not exist in an oblate axisymmetric galaxy. 
We include two SR components, one with the same parameters as 
the NR component, and for the other we set w = 0.5/(— a) and 
u — — 1.0/(— a), and we choose the sense of rotation in the oppo- 
site direction. The latter implies a compact counter-rotating com- 
ponent, which is clearly visible in the kinematic maps shown in the 
top panels of Fig. [14] The inclination is the same value as the po- 
lar angle $ for the triaxial Abel model, i.e. i — 70°, and the mass 
fractions of the three DF components are respectively 20%, 60% 
and 20%. Due to axisymmetry, all maps of the even (odd) veloc- 
ity moments are bi-(anti)-symmetric and the velocity field shows a 
straight zero-velocity curve. The signatures of the counter-rotation 
are similar in the velocity field and /13 (but anti-correlated), and 
result in a decrease of a and an increase of hA in the centre. 



6.4 Recovery of axisymmetric three-integral models 

We now describe the application of our axisymmetric implementa- 
tion of Schwarzschild's method to the observables of the oblate ax- 
isymmetric Abel model of § 16.31 while highlighting the differences 
with the application in triaxial geometry described in Section|5] 



6.4.1 Axisymmetric Schwarzschild model fit to observables of an 
oblate axisymmetric Abel model 

We use the implementation of Schwarzschild's method in axisym- 
metric geometry that is described in detail in Cappellari et al. 
(2006). The main differences with respect to our triaxial implemen- 
tation are certain simplifications due to the extra symmetry. There 
are no twists in the surface brightness and of the four families of 
regular orbits only the short-axis tube orbits are supported. We use 
the same set-up as in the triaxial case, but since there are no box 
orbits, the additional dropping of orbits from the equipotential sur- 
face is not needed. 

Fig. U3] shows that the (simulated) observables of the oblate 
axisymmetric Abel model (top panels) are very well matched by 
the best-fit axisymmetric Schwarzschild model (bottom panels). 
The kinematics of the main body as well as the signatures of the 
counter-rotating core are accurately fitted within the (added) noise. 

6.4.2 Intrinsic velocity moments 

It is convenient to analyse the intrinsic velocity moments of (oblate) 
axisymmetric models in cylindrical coordinates (R, <f), z). Because 
of axisymmetry the models are independent of the azimuthal angle 
4>, and it is sufficient to consider the meridional (R, z)-plane. The 
analysis of the intrinsic velocity moments in the (R, z)-plane is 
similar to that for the triaxial case in the (a;, z)-plane (§ 15. 3b . In this 
case, the mean azimuthal rotation (v^,), perpendicular to the merid- 
ional plane, is the only non-vanishing first-order velocity moment. 
In Fig.[l5] we compare the values of (1)4,) /orms, indicated by the 
colours, for the Abel model (left panel) with those for the fitted 
Schwarzschild model (right panel). The root-mean-square velocity 
dispersion ctrms is defined as orms = { a R + G % + °"z)/3- The 
azimuthal axis of the velocity ellipsoid, with semi-axis length 
defined as = {v^} — {v^,) 2 , is perpendicular to the meridional 
plane. The cross sections with the meridional plane are indicated 
by the ellipses in Fig.[T5] where the semi-axis lengths follow from 
d5.lt by replacing (a;, z) with (R, z). 

As in the triaxial case the density (solid curve) is well fitted by 
the axisymmetric Schwarzschild model (dashed curve). The Abel 
model shows a strong gradient in (v^) / orms, which is correctly 
recovered by the axisymmetric Schwarzschild model. The absolute 
difference is on average less than 0.06, except near the symmetry 
z-axis. This is likely the result of numerical difficulties due to the 
small number of (sampled) short-axis tube orbits that contribute 
in this region. The shape and orientation of the ellipses are nearly 
identical, indicating that the anisotropic velocity distribution of the 
Abel model is reproduced within the expected uncertainties due to 
the errors in the simulated kinematics. The axis ratios ct_/<j+ and 
(j^/<7+ of the velocity ellipsoid are on average recovered within 
~ 5%. 



6.4.3 Three-integral distribution function 

In the oblate axisymmetric case, all (regular) orbits are short- 
axis tube orbits with I2 = \L 2 and energy E ranging from 
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Figure 14. Kinematic maps for an oblate axisymmetric Abel model (top row) and converted to observables with realistic measurement errors added (middle 
row; see § 16. 3t and for the fitted axisymmetric Schwarzschild model (bottom row; see § 16.41 . Parameters and colour scale are as in Fig. [8] 



min [V c fi (A)] to zero. The expression for AV in eq. d5.5t reduces 
to 



AV(E, L„h) = 



(v-A) 



47T 

~\ J J (\ + a)(\+j)(v + a)(v+ 1 ) 

-7 A mi „ 



(X+a)(u+a) 



[E-VaW] [E-V«s(v) 



A\Av (6.9) 



where as before f max , A m in and A max are the solutions of E = 
V e f[{r) (see Fig. 23 of de Zeeuw 1985a). The factor in front of the 
double integral includes the factor 2-n from the integration over the 
azimuthal angle (f>. 

In Fig. [16] we compare the orbital mass weight distribution 
of the input oblate Abel model (top panels), with that of the best- 
fit axisymmetric Schwarzschild model (bottom panels). The three- 



integral mass weight distributions are quite similar, even in the pan- 
els with a relatively low mass content. The average fractional error 
is ~ 30 %, and if we consider in each panel the mass weights above 
the mean value, which together contribute more than half of the 
total mass, the fractional error decreases to around ~ 20 %. Be- 
cause of possible strong point-to-point fluctuations as discussed in 
§ 15.4.31 we also show in Fig. [T7] the orbital mass distribution as 
function of each of the three integrals of motion separately by col- 
lapsing the cube in (E, 12,1a) in the remaining two dimensions. 
Besides the total distribution, we show separately the contributions 
from the NR component and the two opposite rotating SR compo- 
nents in the input oblate Abel model (see § |6.3t . While the com- 
pact counter-rotating SR component (blue) is nicely reproduced by 
the best-fit axisymmetric Schwarzschild model, the mass assigned 
to the main SR component is too high (~ 10 % of the total mass), 
which also results in an underestimation of the NR component. This 
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Figure 15. The mean azimuthal motion (v^) perpendicular to the meridional plane, normalised by (7rms> f° r an oblate axisymmetric Abel model (left) and 
for the best-fit axisymmetric Schwarzschild model (right). Parameters as in Fig. [9] 



V 




■ * 



■ ■ r r J | 



• • I'll Ml 



« * ■ 



Figure 16. The mass weight distribution for the input oblate Abel model (top) and for the best-fit axisymmetric Schwarzschild model (bottom). Parameters as 
in Fig. 1121 The second integral of motion I2 = \L\ > 0, where L z is the component of the angular momentum parallel to the symmetry 2-axis. 



is reflected in the average absolute difference in mass as function 
of E, which is ~ 1.3 %. As for the triaxial case, the recovery for 
I2 and 73 is less good with average uncertainties of ~ 2.1 % and 
~ 2.4 %, respectively. 

A similar good recovery was found by Krajnovic et al. (2005) 



for the case of a two-integral DF f(E,L z ), which implies an 
isotropic velocity distribution in the meridional plane. Thomas et 
al. (2004) showed that their independent axisymmetric numerical 
implementation of Schwarzschild' s method is similarly able to re- 
cover an analytical f(E,L z ). Our results show that the orbital 
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Figure 17. The orbital mass weights for the input oblate Abel model (diamonds connected by solid curves) and for the best-fit axisymmetric Schwarzschild 
model (crosses connected by dotted curves). The parameters are as in in Fig. 1131 except that rotation can only come from short-axis rotating (SR) components, 
for which the two directions of rotation are indicated separately. 



mass weight distribution that follows from a fully three-integral DF 
f(E, L z , I3) can be recovered as well. 



7 DISCUSSION AND CONCLUSIONS 

We have extended the Abel models introduced by DL91 and gener- 
alised by MD99, and shown that, in addition to the intrinsic velocity 
moments, the full LOSVD of these models can be calculated in a 
straightforward way. We have then used the Abel models to con- 
struct realistic axisymmetric and triaxial galaxy models to test the 
accuracy of Schwarzschild' s orbit superposition method. 

Although Abel models have separable potentials with a central 
core and assume a specific functional form for the (three-integral) 
DF, they display a large range of shapes and their observables, 
which can be calculated easily, include many of the features seen 
in the kinematic maps of early-type galaxies. We have used an 
isochrone Stackel potential that in the axisymmetric limit reduces 
to the Kuzmin-Kutuzov model and becomes Henon's isochrone in 
the spherical limit. Because of the simple form of the isochrone 
potential, the resulting Abel models are ideally suited to test nu- 
merical implementations of the Schwarzschild orbit superposition 
method. The calculation of AV, needed when comparing the or- 
bital mass weight distribution of the Schwarzschild models with 
the three-integral DF of the Abel models, simplifies significantly 
for this case. 

Integral-field observations in principle provide the LOSVD 
as a function of position on the sky, so that it is a function 
C{x' , y' , v z i) that depends on three variables. The oblate axisym- 
metric and triaxial galaxy models we have constructed, have a DF 
which is a sum of Abel components f(S) = f(—E + V0I2 + 11I3) 
with different values of the parameters w and u, so that the DF is 
a function of three variables as well, namely the integrals of mo- 
tion E, I2 and J3. We have shown that the simulated integral-field 
observables of these models are matched in detail by the best-fit 
Schwarzschild model. This does not automatically imply that the 
intrinsic velocity moments and the three-integral DF — which are 
not directly fitted — are also correctly recovered. 

First consider three-integral oblate models, i.e., with a DF that 



is a function f(E, L z , I3). In the special case that a galaxy hap- 
pens to be well approximated by a two-integral DF f(E, L z ), the 
density p(R, z) uniquely determines the even part of f(E, L z ) and 
the mean streaming p(v^) in the meridional plane fixes the part 
of f(E,L z ) that is odd in L z (Dejonghe 1986). Ignoring non- 
uniqueness in the deprojection of the surface density E (Rybicki 
1987) and the mean streaming motion V on the plane of the sky, 
these two quantities define a two-integral DF completely. The ob- 
served velocity dispersion and higher-order velocity moments of 
the LOSVD then provide additional information, which for exam- 
ple can be used to constrain the inclination (e.g., Cappellari et al. 
2006). However, the reliability of the derived inclination, of course, 
depends on the correctness of the assumption of a two-integral DF. 
In the more realistic case of a three-integral DF f(E,L z ,Iz), such 
a one-to-one relation with (the velocity moments of) the observed 
LOSVD C(x' ,y' ,v z i) has not been established. Nevertheless, we 
showed that, given integral-field observations of the velocity mo- 
ments of the LOSVD (up to hi), recovery of the full three-integral 
DF is possible with Schwarzschild's method, for the correct incli- 
nation and mass-to-light ratio. 

In the triaxial case, the DF is again a function of three inte- 
grals of motion, but the orbital structure in these models is sub- 
stantially richer than in the oblate axisymmetric models, with four 
major orbit families, instead of only one. This introduces a funda- 
mental non-uniqueness in the recovery of the DF. Whereas in the 
oblate axisymmetric case p(R, z) uniquely defines the even part 
of f(E, L z ), in the (separable) triaxial case the density p(x, y, z) 
does not uniquely determine the even part of f(E, 12,13), although 
both of these are functions of three variables (Hunter & de Zeeuw 
1992). It is (yet) unknown how much specification of C(x' ,y' , v z i) 
can narrow down the range of possible DFs further, even ignor- 
ing the non-uniqueness caused by the required deprojection of the 
surface brightness. Our results show that Schwarzschild's method 
recovers the correct orbital mass weight distribution associated to 
f(E, I2, 13). Given the very large freedom in the orbit choice for 
this case, the modest resolution of our orbit library, and the result- 
ing approximations in the evaluation of the phase space volume, 
the agreement between the orbital mass weights found in § 15.41 is 
in fact remarkable. It may be possible to improve the DF recovery 
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further by refining the sampling of the orbits and the regularisation 
of the orbital mass weights. 

Our analysis shows that it is clear that Abel models are useful 
for testing orbit-based modelling methods such as Schwarzschild's 
method. In particular the oblate limiting case with a Kuzmin- 
Kutuzov potential (§ |6,4t provides a new and convenient test for 
existing axisymmetric Schwarzschild codes. Furthermore, because 
Abel models with a few DF components can already provide quite 
a good representation of observed early-type galaxies, they can be 
used as a way to (numerically) build three-integral dynamical mod- 
els of these galaxies (see e.g. MD96 for an application to Centau- 
rus A). 

We conclude that Schwarzschild's method is able to recover 
the internal dynamical structure of realistic models of early-type 
galaxies. We show in vdB07 that Schwarzschild's method also al- 
lows for an accurate determination of the mass-to-light ratio and 
provides significant constraints on the viewing direction and in- 
trinsic shape. The axisymmetric Schwarzschild method has already 
been successfully applied by us and other groups to determine the 
black hole mass, mass-to-light ratio (profile), dark matter profile as 
well as the (three-integral) DF of early-type galaxies. With our ex- 
tension to triaxial geometry, described in detail in vdB07, we are 
now able to model early-type galaxies — in particular the giant 
ellipticals — which show clear signatures of non-axisymmetry, in- 
cluding isophote twist, kinematic misalignment and kinematic de- 
coupled components. Moreover, since triaxial galaxies may appear 
axisymmetric (or even spherical) in projection, we can investigate 
the effect of intrinsic triaxiality on the measurements of e.g. black 
hole masses based on axisymmetric model fits to observations of 
galaxies. Work along these lines is in progress. 
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APPENDIX A: LIMITING CASES 

When two or all three of the constants a, (5 or 7 that define the 
confocal ellipsoidal coordinate system are equal, the triaxial Abel 
models reduce to limiting cases with more symmetry and thus with 
fewer degrees of freedom. The oblate and prolate axisymmetric 
limits are described in § 16.11 DL9 1 derived the velocity moments 
for the non-rotating Abel models for elliptic discs and in the spher- 
ical limit. We summarise their results and give the rotating Abel 
models as well as the expressions for the LOSVD for these limit- 
ing cases. At the same time, we also derive the properties of the 
non-rotating and rotating Abel models in the limit of large radii. 



confocal elliptic coordinates The relations with {x,y) fol- 

low from those in § l2.1l bv setting z — or equivalently v — —7. 
The two integrals of motion E and I2 are given by 



E 



(vI+vI)+U\\,ij\ 



h = \Ll + Ua - P)vl + (a - /3)x 2 U[X, /x,-a]. 



(Al) 



Al.l Velocity moments 

Choosing the DF as f(E, h) = /(S), with S = -E + w I 2 , the 
velocity moments can be evaluated as 



AtZm(A,/z) 



2i + m+2 
7 l + l, m + l 

rip ii x 



S m ax 



x J T im [S top (A, M ) -S] (i+m)/2 /(S)dS, (A2) 

Smin 



with the terms and h\ defined as 

h T = 1 — (t + a) w, t = A,/i. 



(A3) 



As in the general triaxial case, Smin > Sn m , where the expression 
of the latter is given along the ui-axis (u — 0) in Fig.[T] The acces- 
sible part of the (E, ^-integral space is now a triangle, the top of 
which is S top (A, fi) = —U[X, n] + w(X + a)(n + a)U[X, \i, —a]. 

For the NR components we have that S max = Stop (A, n) and 
T^n = B(!±±, 22±1). Of the two possible orbit families, the box 
orbits have no net rotation and the tube orbits rotate around the 
axis perpendicular to the disc (the z-axis). Since this is similar to 
the short-axis tube orbits in the general triaxial case, we refer to 
the rotating type as the SR type. This SR type reaches the region 
of the accessible integral space (the triangle) for which > at 
jU = ol (or I2 > 0). Therefore, S max = Stop (A, —a) and 



T, S * = 2 



arcsin( v / a^) 

sin' 



)d6, 



where 0,2 is defined as 



02 = 



(A + a)/i M [Stop (A, -a) - S] 



(A4) 



(A5) 



(A - (j,) h ( _ a) [Stop (A, jtt) - S] ' 

The integral dA4t can be evaluated in terms of elementary functions 
(e.g., Gradshteyn & Ryzhik 1994, relations 2.513 on p. 160-162). 
The NR velocity moments /j^(X, fi) vanish when either I or m is 
odd, and the SR velocity moments fif™ (A, fi) vanish when I is odd. 
The latter should be multiplied with (— l) m for net rotation in the 
opposite direction. 

The matrix Q, which transforms the velocity components 
(vx,v^,v u ) to (v x ,v y ,v z ), is that for the prolate case given in 
eq. l |6.4t , but with x = substituted. The sign matrix S, projec- 
tion matrix P and rotation matrix R are the same as for the triaxial 
case given in respectively eqs ( I3,3I h |3~51 >. The polar angle is the 
inclination, $ = i, and the azimuthal angle ip the orientation of the 
infinitesimally thin disc (7 = 0) in the plane z — 0. In the expres- 
sion J3.lt for the misalignment angle ip, the triaxiality parameter 
thus reduces to T = 1 — /3/a, with < f3 < a, bracketing the 
limiting cases of a needle and a circular disc. 



Al Elliptic disc potential 

The two-dimensional analogues of the triaxial Abel models are the 
elliptic Abel discs with Stackel potential Vs(X,n) = U[X,fi] in 



A1.2 Line-of-sight velocity distribution 

Starting with the expression for the stellar (surface) mass density 
E* (x', y') = noo (A, (j,) from eq. JA2t . we derive the Abel LOSVD 
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for the elliptic disc in a similar way as in § 13.3.11 for the triaxial 
cas^]. The cross section of the unit sphere with a plane, reduces 
to the cross section of the unit circle X 2 + Y 2 — 1 with a line 
e\X + e z Y = v z i /g(S). Here, the variable X is defined as 



X = 



hv\ 



with g(S) = h 



2[S top (A, M )-S] 
hxh^ 



(A6) 



and the variable Y follows by interchanging A <-> p. The coeffi- 
cients of the line are ei = \fh\Mz\ /h and 62 = \fh~^M->,2 / h. We 
can write the normalisation h 2 = h\M^ + h^M^, explicitly as 

h = sin i { [1 — (A + a) w] (G cos cp + D sin <p) 2 

+ [1 — (fi + a) w] (Csin ip — D cos tp) 2 } ' 



(A7) 

with G and D defined in eq. i6.5\ . In this way we can write £* 
as a double integral over S and v z i, so that at a given line-of-sight 
velocity the LOSVD becomes 



SupCv) 



f(S) 



•Smin] 



y/2[G(Vj)-S\ 



A(v z ,,S) dS, 



h 



jf(S)A(v z ,,S) sinTydT?, (A8) 




which vanishes when \v z > \ exceeds the 'terminal velocity' v t = 
g(5' m in). The second expression removes the possible singularity 
at the upper limit of S, given by S up = min[G(« z '), Smax], with 



G(v z >) = Stop(A,/i) — hxh^v 2 , /(2h 2 ). 



(A9) 



Hence, ry up is given by sin 2 ry up = [S up - 5 min ]/[G(u z / ) - Smin]- 
Depending on the integral space accessible by the orbits, the value 
of A(v z i , S) is either zero, one or two. 

For the NR components, S max = Stop(A,/n), and, since the 
full integral space is accessible, Anr = 2, independent of S and 
v z i . In the case of a basis function fs(S) as defined in eq. l |2.24b , 
the integral over S can be evaluated explicitly resulting in 



NR 4 S+1 B(S+ 1,6+1) 1 

l~5 = 7= T Lr "z 



Sxnin 



(A10) 



For SR components with maximum streaming, both v 2 > at 
H = —a and u M > 0, which is equivalently to < Y < ^fai, 
where 0,2 is defined in eq. dA5b . The intersection of the above unit 
circle and line provides the following two solutions 



Y± = e 2 [v,,/g(S)] ± eWl-[vz>/g{S)] 2 



(All) 



Given the values of v z i and S, Asr is thus equal to 0, 1 or 2 if for 
respectively none, one or both of the solutions < Y± < ^Jaa. 

The expression for h in eq. JA7| > shows that the LOSVD in 
eq. dA8t is inversely proportional to sin i. For face-on viewing at 
inclination i = 0° the LOSVD reduces to £*(a/, y')6(v z >). Be- 
cause the velocity perpendicular to the disc v z = 0, the face-on 
LOSVD is zero at all line-of-sight velocities, except at v z i = 
when it equals the surface mass density. For edge-on viewing at in- 
clination i — 90°, the LOSVD follows upon substituting y' = in 
eq. JA8t and integrating over the line-of-sight z . For i < 90°, the 



9 Alternatively, one can invert the relations S = Stop ( A, ft) — ^ (h^v 2 + 
hxv^) and v z i = M31 v\ + M32 to find the Jacobian to transform 
from the coordinates {v\,V)i) to (S,v z r). Leaving out the integral over 
v z i yields the same expression for the LOSVD as in eq. )A8t . 



latter integration is not needed, since at each position (x' , y') on the 
plane of the sky there is only a single (unique) point along z' where 
it intersects the infinitesimally thin disc. The edge-on LOSVD and 
also are thus spatially only one-dimensional functions of x', and 
vanish for non-zero i/'-values. 

Further information on elliptic Stackel discs can be found in 
Teuben (1987), de Zeeuw, Hunter & Schwarzschild (1987), and 
Evans & de Zeeuw (1992). 



A2 Large distance limit 

At large radii, A -> r > —a, so that the confocal ellipsoidal 
coordinates of § |2.1| reduce to conical coordinates (r, v), with r 
the usual distance to the origin, i.e., r 2 — x 2 + y 2 + z 2 , and fi 
and v angular coordinates on the sphere. In these coordinates the 
Stackel potential is of the form Vs(r, /x, v) = V(r) + U[fi, v\/r 2 , 
where V(r) is an arbitrary smooth function of r. The corresponding 
integrals of motion are given by 

E = \ {v'i + v y + v 2 z ) + Vs(r, (i, v), 

h = \TL 2 y + \L 2 z + {a-p)^U[^v,-al (A12) 



With the choice J2.9I ) for the DF, the expression for the velocity 
moments becomes 



filmn(r,fi, v) 



2l + m + n+3 



x j T lmn [S top (r,fi,v)- S] {l+m+n+1)/2 f(S)dS, (A13) 

Smin 

where F v and F u are defined as 



(r + a) w — (r + 7) u 

7 — a 



■ v. 



(A14) 



As in the general triaxial case, Smi n > Sn m , where Su m can be 
obtained from Fig. [T] The expressions of S max and T; mn for the 
NR, LR and SR types are those given in ^ I2.3.2if2.3.4l respectivelv. 
but with 5'to P (A, /x, v) (eq. 12. 14t replaced by 



St op (r,n,u) 



-Vs(r,fi,v) 

(u+ a) (v+ a) 
-w — — 

7 — a 

a — 7 



U[fi,v,-a] 

U\fiL,V, -7], 



and the parameters ao and 60 ( 12. 17b reduce to 

S top (r,fi,-f3)-S 
Stop(r,iJ,, v) — S 

(ti + [3)F„ [S top (r,fi,-/3)-S] 



(A15) 



(A 16) 



60 



(H - v) F ( _ /3) [S top (r, H, v ) - S] ' 



which by interchanging v <-> /j, become a\ and 61, and in turn 0,2 
and &2 follow by /3 <-» a. 

In the conversion to observables described in § [3] in the ma- 
trix Q, which transforms the velocity components (« r ,u M ,w y ) to 
(v x ,v y ,v z ), all terms A + a (a = —a, — /3, —7, /i, !/) cancel out 
(cf. eq. 25 of Statler 1994). The expression for the LOSVD follows 
from that of the triaxial case in eq. ( 13.27b by substituting ff M „ = 1, 
H u \ = r 2 F v ,Fxv = r 2 F^ and S top {\, (jl, v) = S top {r,fj,, v). 

Suppose now that at large radii r, the function V(r) in the 
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Stackel potential vanishes and we keep in the above expressions 
only the dominant terms. In this case, F^, F v and Stop reduce to 
functions of /i and v only. As a result, the velocity moments dA 1 3t 
are independent of r , except for the prefactor 1 / r m + n +' 2 ) an d there- 
fore are scale-free. Once we have calculated the velocity moments 
at a radius r, those at radius r' — qr, with q a constant, follow by 
a simple scaling, fj,i mn (r' , fi, v) = /i imn (r, ft, v)/q m+n+2 . The 
same holds true for the line-of-sight velocity moments Hk{r, fj,, v), 
but not for the LOSVD. 

A3 Spherical potential 

When a = /3 = 7, both /1 and v loose their meaning and we 
replace them by the customary polar angle 9 and azimuthal angle <f>. 
The expressions for the Abel models in these spherical coordinates 
(r, 9, (f>) follow in a straightforward way from those in § IA2l for the 
large distance limit in conical coordinates (r, fj,, v). 

The Stackel potential Vs = V(r) is spherically symmetric. 
The expressions for the integrals of motion follow from JA12b , 
where for I2 and Is the right-most terms vanish. The triaxiality 
parameter T is now a free parameter, so that, together with the pa- 
rameters w and u, we can rewrite S = — E + w I2 + u I3 as 

S = -E + \uhl + § [(1 - T)u + Tw]L 2 y + \wL%. (KM) 

This means that with the choice l |2.9t for the DF, we cover the 
most general homogeneous quadratic form in the velocities that is 
allowed by the integrals of motion in a spherical symmetric po- 
tential, i.e., the energy E and all three components of the angular 
momentum vector L (cf. DL91). These include the models consid- 
ered by Osipkov (1979) and Merritt (1985) with the DF of the from 
f(-E ± L 2 /r 2 a ) and those studied by Arnold (1990) with a more 
general DF of the form f(—E ± L 2 /r 2 ± L 2 /r 2 ). These models 
follow by setting u — w — ±2/r 2 , and by taking u — ±2/r 2 , 
w — u ± 2/r 2 and T = 0, respectively. 

A3.1 Velocity moments 

The velocity moments follow from eq. dA13K with 

F T = i- !(«, + «) 

+ i(w - u) [cos 2 9 + T (sin 2 9 sin 2 <j> - 1) ± >/a] , (A18) 
where the positive and negative sign are for F^ and F v , and 

A = [sin 2 9 + T (sin 2 sin 2 - 1)] 2 

+ 4T sin 2 9 cos 2 9 sin 2 <f>. (A19) 

Taking a — f3 = 7 in Fig. Q] we see that the boundaries on w and 
u both vanish. The separatrices L\ and L2, defined in eq. J2.13I ), 
reduce to the negative ui-axis and the line w — u, respectively. 
Furthermore, S max = Stop = ~V(r), and for Ti mn we use the 
expression l |2.15| l. The resulting velocity moments ^j mn (r, 9, tfi), 
which are in general not spherically symmetric, vanish when either 
I, m or n is odd. 

The latter implies no net rotation, which is the case when 
the (conserved) angular momentum vectors L for the orbits are 
randomly oriented. We may introduce net rotation by assuming 
that (a fraction of) the orbits have a preferred sense of rotation 
around an angular momentum vector Lo that points in a specific 
direction given by 9q and <f>o. Using the projection matrix P in 



eq. \3.4\ with $ = 9q and tp — tj>o, we transform to the co- 
ordinate system (r' = r, 9', </)'), in which Lo is aligned with 
the z'-axis. If we next set the DF to zero for L z i < 0, we find 
(ij m „(r, 9', </)') = |/i; mn (r, 0'), which does still vanish when 
I or m is odd, but is non-zero when n is odd, resulting in maxi- 
mum streaming around the z'-axis, and multiplication with (— l) n 
for opposite direction of rotation. With the inverse of the projection 
matrix, we can then transform these velocity moments to the origi- 
nal coordinates system (r, 9, cf>). In this way, we can build spherical 
Abel models, which in addition to a non-rotating part consist of a 
component or several components with a preferred rotation axis. 
Mathieu, Dejonghe & Hui (1996) used this approach to construct a 
dynamical model of Centaurus A, with a spherical potential, but tri- 
axial luminosity density and DF components with rotation around 
the apparent long and short axis. 

From the customary definition of the spherical coordinate sys- 
tem, x = r sin 9 cos <j>,y = r sin 9 sin <j> and z — r cos 9, it follows 
directly that the matrix Q, which transforms the velocity compo- 
nents (v r , v e , ) to (v x , Vy , v z ) , is given by 

(sin 9 cos 4> cos 9 cos 4> — sin <A 
sin 9 sin <j> cos 9 sin (f> cos (f> . (A20) 
cos6> -sin0 / 

In case the orbits have no preferred sense of rotation, we may set 
the viewing angles 1? = ip — without loss of generality, so that 
with ip = from eq. i3.l\ . (x',y') = (y,—x) on the plane of 
the sky and z' — z along the line-of-sight, and similarly for the 
Cartesian velocity components. 



A3. 2 Line-of-sight velocity distribution 

There is no obvious further simplification of the LOSVD for ro- 
tating components. For the NR components, the LOSVD follows 
from eq. J3.27t with S ma x = Stop = — V(r) and A£' = 2tt, 
or from eq. l |3.28t after substituting the basis function fs{S) from 
eq. (12.241 . Since the line-of-sight velocity v z i = sgn(z) [cos 9v r — 
sinOve], it follows that M\ x = cos 2 9, M 32 = sin 2 9 and 
M33 = in eq < f3TT9t for h. Moreover, = 1 and H tX = r 2 F T 
(r = fj,, v), with F T given in eq. JA181 >. 

Carollo et al. (1995) compute the LOSVD for Osipkov-Merritt 
models with f(—E—L 2 /(2r 2 )), which is a special case of the Abel 
DF f(S), that follows from eq. (JaTTJ by setting u = w = — 1/r 2 . 
Substituting the latter in eq. l lA18t , we find that — H\ v — 
(r 2 + r 2 )/r 2 , so that 

h 2 = (r 2 a + r 2 )(r 2 a + r 2 -R' 2 )/r 4 a , (A21) 

gm - -y^)- rl r ;X-n^ (A22) 

with radius R' = r sin 9 on the plane of the sky. After substitution 
in eq. 0.27b . and transforming the integral over Az 1 to dr, we find 
the following LOSVD 

oo Q(v M i) 
C(R',v 2 ,)=4nJ Wr2 r _ Ri2 J /(S)dSdr. (A23) 

R' S min 

This is the same as the (unnormalised) velocity profile in eq. (27) 
of Carollo et al. (1995), with (froo, the lower limit of their (relative) 
potential <E>(r) = —V(r), equal to Smin, which from Fig.[T]in this 
case has Sum = as lower limit. Their function g(r, R') and upper 
limit Qmax in eqs (25) and (26), are equivalent to respectively the 
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Table Bl. The function M for odd s. 



iij M(s,i,j;a,b, 



100 
300 
310 
301 
500 

510 
501 
520 
511 
502 



|(4- a- 6) + \(b - a) sin 20 
— \<j>— j sin 20 
-\<t> + | sin 20 

| (24 - 12a - 126 + 3a 2 + 36 2 + 2ab) <j> 

+|(6 - a)(3 - a - 6) sin 20 + i(6 - a) 2 sin 40 

-|(6 - 3a - 6)0- |(3 - 2a) sin 20 - ±(b - a) sin 40 

-|(6- 36-a)0 + |(3-26) sin 20+ ^(b-a) sin 40 

|0+ ± sin 20 + ^ sin 40 

j4>- jg sin 40 

|0- | sin 20+ i sin 40 



inverse of h and G(v z /) in eqs JA2U and d A22b above. The well- 
known isotropic case follows upon taken the limit r a — ► oo, so that 

f(S) -» /(-£), ft - 1 and G(tv) - -V(r) - v 2 z ,/2. 



APPENDIX B: THE FUNCTION M 

The function M that appears in the velocity moments of the rotat- 
ing Abel components is defined as 



4> 



. . . . aV/ay i-^[i-P(9)r +1 Jfl 

M(s, h r,a,b^) = J (^J — dfl, 



(Bl) 

with = a cos 2 8 + 6 sin 2 9. For odd s, corresponding to odd 
velocity moments, the integral can be evaluated in a straightforward 
way in terms of elementary functions. In Table IB 1 1 we give the 
resulting expressions for s = 1, 3, 5. 

For even s, the integral can be evaluated in terms of the (in- 
complete) elliptic integrals. To simplify the numerical evaluation 
we use Carlson's (1977) symmetrical forms Rf, Rd and Rj (for 
the relations between both forms see e.g. de Zeeuw & Pfenniger 
1988). In Table lB2l we give the expressions for s — 0, 2, 4, where 
we have introduced the following quantities based on these sym- 
metric elliptic integrals 



F 
D 

J 



VI - a, sin a, a1 
itf+cos (f>, A ,1) 



-Rd(cos 2 0, A 2 , 1), 



3%/l-a 

(b - a) sin 3 . 2 , A 2 , 

— , 7= -Rj[cos 4>,A ,1, — 

3a 2 vl — a a 



(B2) 



with A 2 = [1 — p(0)]/(l — a), and we have defined the terms 



A 



P = 




i / b 

arctan \ \ — tan< 
^06 \ V a 



sm cos 



(B3) 



In Fig. IB II we show the M (s, i , j; a, b, 0) as function of <j> for the 
case that a = 0.1 and b — 0.5, up to order s = 5. 

We now consider some special cases. When either a or b is 
zero, the corresponding velocity moments vanish (eqs 12.161 and 



12. 19t , and when Oj > bj the arguments of the function .M are 
interchanged (eqs 12. 18112.201 and 12.2 It . This means we only have 
to consider the range < a < b, together with < <j> < tt/2, 
since A4 vanishes when = 0. 

When a = b, it follows that p(6) = a in eq. (Bll , so we can 
separate M(s, a, a, <j>) — Mi(s, a) A42(i,j; <j>), where 



Mi(s,i,j;a) = 
M 2 (i,j;(t>) = 



(B4) 



o 



sin 2j 060. 



For a = 1, the expression for A4i simplifies to ( — l) 1+ -' (j + j)!. 
The integral in the expression for M2 can be evaluated explicitly 
using e.g. the relations 2.513 of Gradshteyn & Ryzhik (1994). For 
tf> = tt/2, it reduces to the beta function B(i + 1/2, j + 1/2). 

When o < 6 = 1, the elliptic integrals become elementary, so 
that the quantities F, D and J in eq. dB2t reduce to 



, 7T 

tan - + - 
1 4 2 



D = 



a sin 

■ r — 



I — a y/l — a 



(B5) 



J = F 



arctan 



1 - a 



Although F diverges when — > 7r/2, substitution of these reduced 
quantities in the expressions of M for even s (Table lB2l . shows 
that all terms with F cancel. For = 7r/2, the function A4 is thus 
everywhere finite, with A — tt/(2Vot) and P = Q = 0. 



APPENDIX C: EDGEWORTH EXPANSION 

For the (re)construction of the LOSVD from its true line-of-sight 
velocity moments, one can use the well-known Gram-Charlier se- 
ries, the terms of which are simple functions of the true moments 
(see e.g. Appendix B2 of van der Marel & Franx 1993), but it has 
poor convergence properties. The terms in the Edgeworth (1905) 
expansion are also directly related to the true moments, but since 
it is a true asymptotic expansion its accuracy is controlled, so that, 
unlike the Gauss-Hermite and Gram-Charlier expansions, conver- 
gence plays no role (see Blinnikov & Moessner 1998 for a compar- 
ison between the expansions and for further references). 

The Edgeworth expansion of the LOSVD up to order N is 
given by 



27TCT 



1 +Y. Dn 



with w = (v — V)/(T and 

- 1 

Du= W " + 2(i-l)Hn 1 

0,-2} <=3 2-2 



(CI) 



(C2) 



The Hermite polynomials TL m are related to those defined by van 
der Marel & Franx (1993) as H m {w) = \fm\ H m (w/^/2). We 
have defined / = Yl^Zi lj> wnere me sel;s {h } are me non-negative 
integer solutions of the Diophantine equation 

lj + 2lj + h in - 2)/ n _ 2 = n - 2, n > 3, (C3) 
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Table B2. The function JV[ for even s. 



iij A4(s,i,j;a,b, 
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Figure Bl. The function M{x. a, b, <j>) defined in eq. IB It plotted against <f>, for a = 0.1 and b = 0.5, up to order s = 5. The curves in the left two 
panels are for odd values of s corresponding to the odd velocity moments, whereas the curves in the right two panels are for even values of s. The indices of 
the labels M ai j refer to the first three parameters of the function A4. 



Substituting these solutions, we find up to order N — 5 



to the true moments respectively as 



/.ED, n vi e 



2iva 



(C4) 



The lower-order moments E, V and a are equivalent to those in 
eq. 0.4U . while the higher-order moments di (i > 3) are cumu- 
lants of the true moments 



di = 



{i k } k=i 



(C5) 



so that 

da =6, d 4 =6-3, and d B =&-10£i. (C6) 
The central moments £1 (skewness), £2 (kurtosis) and £3 are related 



(Moc) £l = Mo M3 - 3/Lto /X1/U2 + 2nx, 

(Hoa) 4 & = /io/i4 — 4/io/ii ^3 + 6^0 ^i?/i2 — 3/ii, 



(Ho<t) S £3 = Mo Ms - 5 /Lto Mi M4 + 10 Mo Mi 
— 10 Mo M? M2 + 4 Mi- 



M3 



(C7) 
(C8) 

(C9) 



Substituting the line-of-sight true moments Mfc f° r fe = 0, . . . ,K, 
we can compute £™ (w) at each position on the plane of the sky. 

In Fig. |C1| we show an example of a LOSVD (black solid 
line) computed directly via eq. J3.27t for the triaxial Abel model 
constructed in § 14.31 The Edgeworth LOSVD (red solid line) is 
constructed from the true line-of-sight velocity moments, based 
on the intrinsic velocity moments computed via eq. d2. 10t . The 
Edgeworth reconstruction approximates (very) well the directly- 
computed LOSVD, as well as the corresponding best-fit Gauss- 
Hermite series (blue dashed line). In Fig. IC2I we show the Gauss- 
Hermite moments after fitting at each (aperture) position on the 
sky-plane the directly-computed LOSVD (top panels) as well as the 
reconstructed Edgeworth LOSVD (bottom panels). The resulting 
maps are very similar, except for a suppression of the higher-order 
Gauss-Hermite moments in case of the Edgeworth reconstruction. 
This is expected, since intrinsic velocity moments of order higher 
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Figure CI. Line-of-sight velocity distribution (LOSVD) of the triaxial Abel 
model presented in § 14.31 at the (aperture) position on the sky-plane given at 
the top of the figure. The black solid curve is the LOSVD computed directly 
via eq. B.27I . The red solid curve show the Edgeworth LOSVD constructed 
from the true line-of-sight velocity moments, based on the intrinsic velocity 
moments computed via eq. j2.10t . The Gaussian and the higher order terms 
of the Edgeworth expansion 4C H are shown by the red dotted curves. The 
blue dashed curve shows the best-fit Gauss-Hermite LOSVD. 



than N = 5 are needed to accurately determine the wings of the 
LOSVD. Nevertheless, this comparison is important to show the 
correctness of both (independent) approaches, and that the Edge- 
worth expansion provides a reliable and efficient way to reconstruct 
the LOSVD (and obtain Gauss-Hermite moments) from true mo- 
ments that are in general (numerically) easier to compute than the 
full LOSVD. 

This paper has been typeset from a Tj^X/ LTgX file prepared by the 
author. 
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Figure C2. Maps of the surface mass density (£; in 10 4 M pc~ 2 ), mean line-of-sight velocity V and dispersion a (both in km s 1 ), and higher order 
Gauss-Hermite moments /13 and /14, of the triaxial Abel model constructed in § 14.31 The top panels follow from fitting, at each (aperture) position on the 
sky-plane, Gauss-Hermite series to the LOSVD computed directly via eq. )3.27> . For the bottom panels the fit is applied to the Edgeworth LOSVD constructed 
from the true line-of-sight velocity moments, based on the intrinsic velocity moments computed via eq. 12.10) . 



